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Abstract 

A  central  characteristic  feature  of  an  important  class  of  hyperbolic  PDEs  in  odd- 
dinrension  spaces  is  the  presence  of  lacunae ,  i.e.,  sharp  aft  fronts  of  the  waves,  in 
their  solutions.  This  feature,  which  is  often  associated  with  the  Huygens’  principle, 
is  employed  to  construct  accurate  non-local  artificial  boundary  conditions  (ABCs) 
for  the  Maxwell  equations.  The  setup  includes  the  propagation  of  electromagnetic 
waves  from  moving  sources  over  unbounded  domains.  For  the  purpose  of  obtaining 
a  finite  numerical  approximation  the  domain  is  truncated,  and  the  ABCs  guarantee 
that  the  outer  boundary  be  completely  transparent  for  all  the  outgoing  waves.  The 
lacunae-based  approach  has  earlier  been  used  for  the  scalar  wave  equation,  as  well 
as  for  acoustics.  In  the  current  paper,  we  emphasize  the  key  distinctions  between 
those  previously  studied  models  and  the  Maxwell  equations  of  electrodynamics,  as 
they  manifest  themselves  in  the  context  of  lacunae-based  integration.  The  extent  of 
temporal  nonlocality  of  the  proposed  ABCs  is  fixed  and  limited,  and  this  is  not  a 
result  of  any  approximation,  it  is  rather  an  immediate  implication  of  the  existence  of 
lacunae.  The  ABCs  can  be  applied  to  any  numerical  scheme  that  is  used  to  integrate 
the  Maxwell  equations.  They  do  not  require  any  geometric  adaptation,  and  they 
guarantee  that  the  accuracy  of  the  boundary  treatment  will  not  deteriorate  even  on 
long  time  intervals.  The  paper  presents  a  number  of  numerical  demonstrations  that 
corroborate  the  theoretical  design  features  of  the  boundary  conditions. 
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1  Introduction 


Sharp  aft  fronts  of  the  waves,  or  lacunae,  is  a  unique  feature  of  certain  hy¬ 
perbolic  solutions  in  odd-dimension  spaces.  The  corresponding  class  of  PDEs 
includes  the  Maxwell  equations  of  electrodynamics  studied  in  the  current  pa¬ 
per,  along  with  the  acoustics  system  (linearized  Euler  equations)  and  the  scalar 
wave  equation,  that  have  been  analyzed  previously,  see  [1-3] .  The  key  idea  of 
using  lacunae  for  computations  is  very  simple.  If  the  sources  of  waves  are  com¬ 
pactly  supported  in  space  and  time,  and  if  the  domain  of  interest  has  finite 
size,  then  it  will  completely  fall  inside  the  lacuna  once  a  certain  time  interval 
has  elapsed  since  the  inception  of  the  sources.  This  interval  can  be  easily  evalu¬ 
ated  ahead  of  time,  and  for  all  subsequent  moments  of  time  the  solution  on  the 
domain  of  interest  will  be  zero,  as  this  domain  will  already  be  behind  the  aft 
front  of  the  waves.  Therefore,  no  computation  will  need  to  be  run  any  further, 
which  puts  a  cap  on  the  long-term  error  build-up.  Moreover,  during  this  prede¬ 
termined  finite  interval  no  wave  can  travel  more  than  a  certain  distance  away 
in  space  from  the  domain  of  interest.  Consequently,  a  finite  auxiliary  com¬ 
putational  domain  can  be  set  up,  which  facilitates  construction  of  artificial 
boundary  conditions  (ABCs).  In  general,  the  term  “ABCs”  is  associated  with 
a  fairly  large  group  of  methods  employed  for  solving  infinite-domain  problems 
on  the  computer.  For  the  purpose  of  obtaining  a  finite  numerical  approxima¬ 
tion  the  original  domain  is  truncated,  and  the  ABCs  provide  a  required  closure 
for  the  resulting  truncated  formulation.  The  literature  on  the  subject  of  ABCs 
is  broad,  and  we  refer  the  reader  to  the  review  papers  [4-6]. 

It  is  generally  known  that  highly  accurate  ABCs  are  nonlocal.  In  particular, 
exact  ABCs  are  always  nonlocal  in  multidimensional  settings.  For  unsteady 
problems,  this  also  implies  nonlocality  in  time.  Nonlocal  ABCs  may  incur 
high  computational  costs  and  require  elaborate  implementation  strategies.  As 
such,  local  alternatives  obtained  either  independently  or  as  an  approximation 
of  nonlocal  ABCs  conditions,  becomes  viable.  These  ABCs  are  usually  inex¬ 
pensive  and  easy  to  implement,  but  may  often  lack  computational  accuracy. 

The  ABCs  for  the  Maxwell  equations  that  we  propose  hereafter  are  nonlo¬ 
cal,  and  their  accuracy  can  always  match  that  of  the  interior  discretization. 
However,  the  extent  of  their  temporal  nonlocality  is  limited,  and  does  not 
increase  as  the  time  elapses.  The  bound  on  temporal  nonlocality  comes  as  a 
consequence  of  the  presence  of  lacunae,  i.e.,  sharp  aft  fronts  of  the  waves,  in 
the  corresponding  solutions.  In  the  literature,  existence  of  the  lacunae  is  often 
referred  to  in  connection  with  the  Huygens’  principle. 

The  lacunae-based  ABCs  are  not  limited  to  the  aforementioned  case  of  the 
compactly  supported  sources  of  waves.  The  case  of  continuously  operating 
sources,  or  in  general,  that  of  a  certain  process/phenomenon  confined  to  a 
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bounded  region  that  manifests  itself  by  radiation  of  electromagnetic  waves  in 
the  far  field,  is  treated  by  employing  special  partitions  along  the  time  axis. 

Several  other  nonlocal  ABCs’  methodologies  for  unsteady  waves  have  recently 
been  put  forward,  most  notably  [7-14];  see  also  the  survey  [15].  In  comparison, 
our  approach  has  a  number  of  distinctive  characteristics,  besides  the  restricted 
temporal  nonlocality  that  does  not  come  at  a  cost  of  accuracy.  The  lacunae- 
based  ABCs  are  obtained  directly  for  the  discretization,  and  can  supplement 
any  consistent  and  stable  scheme.  In  other  words,  they  bypass  the  common 
stage  of  deriving  the  continuous  ABCs  and  subsequently  approximating  them, 
see  [5].  The  lacunae-based  ABCs  withstand  long-term  numerical  integration 
with  no  deterioration  of  accuracy  and  are  not  restricted  geometrically  to  any 
particular  shape  of  the  external  artificial  boundary.  They  also  allow  one  to 
analyze  a  variety  of  cases,  including  the  one  of  moving  sources  of  waves. 

Even  though  the  current  ABCs  for  electromagnetics  are  related  to  those  devel¬ 
oped  previously  for  the  scalar  wave  equation  [1,2]  and  for  acoustics  [3],  there 
are  not  only  similarities  but  also  very  important  differences  between  the  re¬ 
spective  models.  Foremost,  the  source  terms  that  drive  the  Maxwell  equations 
are  not  independent,  they  are  related  via  the  continuity  equation.  The  latter 
must  be  enforced  on  all  stages  of  building  the  lacunae-based  algorithm.  In 
general,  the  current  ABCs’  approach  fits  into  the  universal  framework  of  [16]. 

The  scope  of  numerical  demonstrations  that  we  present  is  limited  to  the  trans¬ 
verse  magnetic  (TM)  mode  with  cylindrical  symmetry.  This  allows  us  to  take 
advantage  of  the  three-dimensional  effects  in  a  two-dimensional  setting. 


2  Lacunae  of  the  Wave  Equation 


Consider  a  three-dimensional  wave  equation,  x  =  (xi1X2,xf)\ 


1  d2cp 
c 2  dt2 


d2p>  d2p>  d\p\ 
dx2  dx2  dx2 ) 


f(x,t),  t>  0 


<P 


t= o 


0 


(la) 

(lb) 


with  homogeneous  initial  conditions.  For  every  ( x,t ),  the  solution  p  =  p(x,t) 
of  problem  (la),  (lb)  is  given  by  the  Kirchhoff  integral  (see,  e.g.,  [17]): 


ip(x,t) 


illl 

Q<Ct 


(2) 
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where  £  =  (6,6,6),  Q  =  I*  -  £|  =  \/ (afi  -  6)2  +  (^2  -  6)2  +  (^3  -  6)2, 
and  d$,  =  dfid^df-j.  If  we  assume  that  the  right-hand  side  (RHS)  f(x,t ) 
of  equation  (la)  is  compactly  supported  in  space-time  on  the  domain  Q  C 
M3  x  [0,  +00),  then  formula  (2)  immediately  implies  that 


0  for  (x,t)  G  P)  {(&,£) 

(S.ReQ 


x  —  £|  <  c(f 


r),  t  >  r} 


(3) 


The  region  of  space-time  defined  by  formula  (3)  is  called  lacuna  of  the  so¬ 
lution  Lp  =  (p(x,t).  This  region  is  obviously  obtained  as  the  intersection  of 
all  characteristic  cones  of  equation  (la)  once  the  vertex  of  the  cone  sweeps 
the  support  of  the  RHS:  supp  /  C  Q.  From  the  standpoint  of  physics,  lacuna 
corresponds  to  that  part  of  space-time,  on  which  the  waves  generated  by  the 
sources  f(x,t),  supp  /  C  Q,  have  already  passed,  and  the  solution  has  be¬ 
come  zero  again.  The  phenomenon  of  lacunae  is  inherently  three-dimensional. 
The  surface  of  the  lacuna  represents  the  trajectory  of  aft  (trailing)  fronts  of 
the  waves.  The  existence  of  aft  fronts  in  odd-dimension  spaces  is  known  as 
the  Huygens’  principle,  as  opposed  to  the  so-called  wave  diffusion  which  takes 
place  in  even-dimension  spaces,  see,  e.g.,  [17].  Note  that  the  aft  fronts  and  the 
lacunae  would  still  be  present  in  the  solution  of  the  wave  equation  (la)  if  the 
homogeneous  initial  conditions  (lb)  were  replaced  by  some  inhomogeneous 
initial  conditions  with  compact  support. 

The  notion  of  lacunae  (or  lacunas)  was  first  introduced  and  studied  by 
Petrowsky  in  [18]  (see  also  an  account  in  [19,  Chapter  VI])  for  a  variety  of 
hyperbolic  equations  and  systems;  and  general  characterization  of  their  coef¬ 
ficients  was  provided  that  would  guarantee  existence  of  the  lacunae.  However, 
since  work  [18]  no  constructive  examples  of  either  equations  or  systems  have 
been  obtained,  for  which  lacunae  would  be  present  in  the  solutions,  besides 
the  actual  wave  equation  (la),  as  well  as  those  equations  that  either  reduce 
to  or  are  derived  from,  the  wave  equation. 

The  group  of  numerical  algorithms  for  integrating  the  Maxwell  system  of  equa¬ 
tions  that  we  describe  hereafter  will  be  essentially  based  on  the  presence  of 
lacunae.  However,  the  Kirchhoff  formula  (2)  will  never  be  used  as  an  actual 
part  of  the  algorithm,  it  will  only  be  needed  at  the  theoretical  stage,  for  deter¬ 
mining  the  shape  of  the  lacunae  that  will  later  be  incorporated  into  a  purely 
finite-difference  context.  Previously,  the  idea  of  using  the  Huygens’  principle 
for  constructing  the  ABCs  was  promoted  by  Ting  &  Miksis  [20]  and  Givoli  & 
Cohen  [21].  Both  papers,  however,  have  suggested  to  use  numerical  quadra¬ 
tures  to  approximate  the  integral  (2),  and  then  couple  it  with  the  interior 
solver.  Moreover,  the  approach  of  [20]  has  never  been  actually  implemented 
in  a  practical  computational  setting,  whereas  the  approach  of  [21]  required 
artificial  dissipation  to  be  added  to  the  scheme  to  fix  the  arising  instabilities. 
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3  The  Maxwell  System  of  Equations 


The  propagation  of  electromagnetic  waves  in  vacuum  is  governed  by  the 
Maxwell  equations: 


ldH 

c  dt 
Id  E 

c  dt 


+  V  x  E 
-  V  x  H 


V  •  H  =  0 

V  •  E  =  47rp 


(4) 


driven  by  the  electric  charges  p  =  p(x,t)  and  currents  j  =  j(x,t).  In  system 
(4),  E  is  the  electric  field,  H  is  the  magnetic  field,  c  is  the  speed  of  light,  and 
the  normalization  is  chosen  so  that  both  the  permittivity  and  permeability 
of  vacuum  are  equal  to  one:  £q  —  p,o  —  1.  We  note  that  unlike  the  acoustics 
system  studied  in  [3],  all  the  unknown  quantities  in  (4)  are  vectors. 


System  (4)  contains  a  pair  of  unsteady  equations  and  a  pair  of  steady-state 
equations.  Taking  the  operation  curl  =  Vx  of  one  equation  from  the  first 
pair,  differentiating  the  remaining  unsteady  equation  with  respect  to  time, 
substituting  into  one  another,  and  using  the  identity  2  V  x  V  x  V  =  —  A  V  + 
V(V  •  V)  along  with  the  corresponding  steady-state  equation,  we  arrive  at 
the  following  individual  equations  for  the  fields: 


1  d2H 
c2  dt 2 
1  d2E 
c2  dt 2 


-AH 


-A  E 


4vr 

— V  x 
c 


3 


—47 r 


dj_ 
c2  dt 


+  Vp 


(5) 


Each  of  the  two  equations  (5)  is  a  vector  wave  equation.  Therefore,  if  the 
currents  and  charges  in  system  (4)  were  compactly  supported  in  space  and 
time,  as  in  Section  2,  then  the  RHSs  in  equations  (5)  would  be  compactly 
supported  as  well,  and  consequently,  the  solutions  E  =  E(x,t )  and  H  = 
H  (x,t)  would  have  lacunae  of  the  same  shape  as  prescribed  by  formula  (3). 


At  this  stage,  another  important  difference  compared  to  the  acoustics  system 
of  equations  needs  to  be  outlined.  In  acoustics  [3],  we  have  required  that  the 
velocity  field  be  conservative.  It  was  necessary  and  sufficient  for  the  individual 
equation  for  acoustic  velocity  to  become  a  wave  equation.  ^  In  its  own  turn, 
the  fact  that  a  given  quantity  is  governed  by  the  wave  equation  is  sufficient 
for  the  solutions  driven  by  compactly  supported  sources  to  have  lacunae. 


In  contradistinction  to  that,  individual  equations  that  govern  the  electric  and 
magnetic  fields  always  reduce  to  the  wave  equation,  see  (5).  As  such,  existence 

2  V  can  be  any  vector. 

3  The  other  acoustic  variable,  pressure,  always  satisfies  the  wave  equation,  [3] . 
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of  the  lacunae  for  compactly  supported  sources  is  always  guaranteed.  This, 
however,  comes  “at  a  price.’’  Indeed,  there  are  two  types  of  sources  in  acoustics, 
the  volume  velocity  (scalar)  that  drives  the  continuity  equation,  and  the  force 
(vector)  that  drives  the  momentum  equation.  These  two  sources  are  completely 
independent,  see  [3]  and  also  [22],  In  electrodynamics,  this  is  not  the  case. 
The  currents  and  charges  that  drive  the  Maxwell  equations  (4)  are  always 
connected  via  the  continuity  equation  [23]: 

!  +  VO=0  (6) 

Equation  (6)  may  present  a  certain  limitation  when  building  the  ABCs  for  the 
Maxwell  equations  (4).  Indeed,  our  general  approach  [1,2]  involves  construc¬ 
tion  of  the  special  auxiliary  sources  that  have  to  be  compactly  supported  in 
both  space  and  time,  i.e. ,  in  the  sense  of  Section  2.  It  is  clear,  however,  that 
specifying  a  compactly  supported  current  j  (x,  t )  may  yield  the  distribution  of 
charges  p(x,t)  that  is  not  necessarily  compactly  supported  in  time  because  of 
equation  (6).  For  the  numerical  examples  of  Section  8  this  limitation  does  not 
manifest  itself  clue  to  the  cylindrically  symmetric  setup  that  we  have  adopted, 
see  Sections  5  and  6.  In  general,  however,  additional  attention  may  be  required 
toward  this  issue.  If,  on  the  other  hand,  we  knew  ahead  of  time  that  the  orig¬ 
inal  currents,  as  well  as  charges  in  (4),  were  compactly  supported,  then  the 
RHSs  of  (5)  would  be  compactly  supported  as  well,  i.e.,  equation  (6)  would 
carry  no  restriction  from  the  standpoint  of  existence  of  the  lacunae. 


4  Decomposition  of  the  Original  Problem 


As  indicated  in  Section  1,  artificial  boundary  conditions  are  needed  to  truncate 
the  problem  which  is  originally  formulated  on  an  unbounded  domain,  and  as 
such  enable  its  solution  on  the  computer.  Suppose  that  the  original  formulation 
of  the  problem  involves  the  entire  space  M3,  but  we  are  only  interested  in 
finding  a  fragment  of  the  overall  solution  defined  on  a  given  domain  S(t). 
This  domain  has  a  fixed  shape  and  finite  diameter,  diarn  S(t)  =  d  <  oo,  for 
all  t  >  0,  but  is  allowed  to  move  across  M3  according  to  the  general  law: 

u0  =  Uo(t),  Xq  =  x o(t)  =  I  u0(r)dT  (7) 

Jo 

where  u0  and  x$  are  the  velocity  vector  and  coordinates  of  a  given  point  inside 
S(t),  respectively.  The  only  natural  limitation  is  that  max(  \uo(t)  \  =  k  <  c. 

While  not  making  any  specific  assumptions  regarding  the  nature  of  the  phe¬ 
nomena/processes  that  are  going  on  inside  Sit),  we  assume  that  outside  Sit), 
i.e.,  Wt  >  0  and  x  qL  S(t),  the  appropriate  model  would  be  based  on  the 
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homogeneous  Maxwell  equations  [cf.  equations  (4)]: 

-?  +  Vx£  =  0  V-H  =  0 

idE  (§) 

x  H  =  0  V-E  =  0 

c  dt 

We,  of  course,  assume  that  the  overall  problem,  i.e.,  the  interior  one  that  we 
do  not  elaborate  on,  combined  with  the  exterior  one,  which  is  governed  by 
system  (8),  is  uniquely  solvable  on  M3.  In  other  words,  our  model  may  include 
some  possibly  complex  phenomena  confined  to  the  bounded  region  S(t)  that 
manifest  themselves  by  the  radiation  of  electromagnetic  waves  in  the  far  field. 
The  objective  is  to  actually  solve  the  problem  on  the  domain  S(t)  using  a 
numerical  method,  but  truncate  all  of  its  exterior  replacing  it  with  the  ABCs 
on  the  external  boundary  dS(t).  The  ideal  or  exact  ABCs  would  make  the 
foregoing  replacement  equivalent,  which  means  that  the  solution  obtained  this 
way  would  coincide  on  S(t)  with  the  corresponding  fragment  of  the  original 
infinite-domain  solution.  The  latter,  of  course,  is  not  actually  available  because 
the  far-field  propagation  of  electromagnetic  waves  governed  by  (8)  cannot  be 
computed  directly  at  an  acceptable  cost. 

Let  us  emphasize  that  from  the  viewpoint  of  an  observer  inside  S(t),  the  role  of 
the  ABCs  at  dS(t)  is  only  to  guarantee  that  this  boundary  behave  exactly  as 
if  the  domain  S(t)  were  surrounded  by  vacuum.  In  particular,  the  boundary 
dS(t)  may  not  reflect,  fully  or  partially,  any  outgoing  waves.  Therefore,  as 
far  as  the  ABCs  are  concerned,  we  indeed  do  not  have  to  be  very  specific 
regarding  the  nature  of  the  problem  inside  S(t),  provided,  of  course,  that  the 
overall  interior /exterior  problem  does  have  a  unique  solution. 

The  latter  assumption  is  obviously  of  key  importance.  However,  justifying  it  for 
every  particular  setting  is  beyond  the  scope  of  the  current  paper.  A  relatively 
easy  case  would  be  that  of  linear  scattering  from  a  known  (sufficiently  simple) 
shape  with  known  material  characteristics.  Otherwise,  the  model  inside  S(t) 
may  be  more  sophisticated,  possibly  nonlinear.  In  any  event,  for  the  purpose 
of  constructing  the  ABCs,  the  overall  solvability  will  hereafter  be  assumed, 
including  the  case  when  S(t)  is  moving.  We  only  mention  that  in  the  context 
of  electromagnetics  the  capability  of  the  method  to  handle  moving  sources  of 
waves  is  probably  less  important  than  in  acoustics,  see  [3],  for  the  obvious 
reason  that  all  feasible  speeds  will  be  much  slower  than  the  speed  of  light.  We, 
however,  still  include  this  capability  in  the  design.  One  reason  is  that  there 
may  be  future  applications,  for  which  it  will  be  important,  e.g.,  in  astrophysics. 
Another  reason  is  that  in  electromagnetics,  the  structure  of  the  field  does 
depend  on  the  frame  of  reference,  in  which  this  field  is  analyzed,  see  [23].  Again, 
one  can,  in  principle,  think  of  situations,  in  which  employing  a  moving  frame 
of  reference,  which,  in  turn,  may  imply  moving  sources,  could  be  beneficial 
because  of  a  potentially  simpler  structure  of  the  associated  field. 
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The  first  stage  in  constructing  the  ABCs  will  be  decomposition  of  the  original 
infinite-domain  problem  into  two  sub-problems.  An  auxiliary  problem  will  be 
formulated  that  will  have  the  exact  same  solution  outside  S(t)  as  the  original 
problem  does,  but  will  be  linear  throughout  the  entire  space  M3  and  will  be 
driven  by  the  special  source  terms  concentrated  only  inside  S(t).  Next,  the 
auxiliary  problem  will  be  solved  using  the  lacnnae-based  methodology  [1-3] 
adopted  here  for  the  Maxwell  equations,  and  its  solution  obtained  right  outside 
S(t)  will  be  used  to  supply  the  boundary  conditions  for  the  interior  problem 
solved  inside  S(t).  In  practice,  the  two  aforementioned  stages  will  be  meshed 
together  so  that  both  the  interior  problem  and  the  auxiliary  problem  are 
time-marched  concurrently.  In  so  doing,  the  interior  solution  will  be  used  to 
generate  sources  for  the  auxiliary  problem,  and  the  auxiliary  solution  will  be 
used  to  provide  the  missing  boundary  data  for  the  interior  problem.  The  entire 
algorithm  will  be  implemented  directly  on  the  discrete  level. 


Let  us  define  a  smaller  subdomain  S£(t)  C  S(t)  such  that  x  e  S£(t)  iff  x  e  S(t) 
and  dist(cc,  dS(t ))  >  e,  where  £  >  0.  Let  us  also  introduce  a  multiplier  function 
H  =  fi(x,t)  that  is  smooth  across  the  entire  space  and  Vt  >  0  satisfies: 


f  0,  x  G  Se(t) 

M*,*)  =  |l,  x  &  S(t)  (9) 

U(°,l)>  xeS(t)\Se{t) 


The  curvilinear  strip  S(t)\S£(t)  of  width  e  adjacent  to  the  boundary  dS(t) 
of  the  domain  S(t)  from  inside  will  hereafter  be  called  the  transition  region. 
Assume  that  the  solution  to  the  original  problem  E(x,t ),  H(x,t )  is  known 
on  M3  for  t  >  0,  and  that  E(x,  0)  =  0  and  H(x,  0)  =  0  for  x  <0  Se( 0),  which 
means  that  nothing  is  going  on  outside  the  domain  of  interest  before  t  —  0. 
We  multiply  this  solution  on  the  entire  space  by  p(x,t)  of  (9)  and  obtain: 

E{x,t)  i — *  fji{x,t)E(x,t)  =  E(x,t) 

H(x,t)  i — ■>  fj,(x,i)H(x,t)  =  H(x,t ) 


We  emphasize  that  in  formulae  (10)  we  do  not,  in  fact,  need  to  know  E(x,t) 
and  H(x,t)  on  S£(t),  because  the  multiplier  fi  is  equal  to  zero  there  anyway. 
Moreover,  multiplication  (10)  will  not  change  any  quantities  on  M.3\S(t).  With 
that  in  mind,  we  apply  the  differential  operators  from  the  left-hand  sides  of 
system  (8)  to  the  modified  solution  E(x,t),  H(x,t),  see  (10),  and  obtain: 


in  -  1  dH  -  .  „  Y7  ZT 

- Jm:=  --  +  V  x  E  4npm  :=  V  •  iT 

c  c  at 

-—j  :=~-Vxff  4v rp:=V-E 

c  c  at 


(ii) 


Relations  (11)  define  the  auxiliary  sources  jm(x,t),  j(x,t),  pm(x,t),  and 


p{x,t)  on  the  space  M3  for  t  >  0.  Clearly,  for  x  G  Se(t)  and  t  >  0  we  have: 

jm{x,  t)  =  0  pm(X,t)=  0 

j(x,t)  =  0  p(*,t)  = 0 

because  p(x,t)  =  0  for  a;  e  S^f).  Equalities  (12)  also  hold  for  cc  e  M3\5'(f) 
and  all  t  >  0,  because  the  modified  quantities  E(x,t)  and  H(x,t)  of  (10)  co¬ 
incide  with  E(x,t)  and  H(x,t),  respectively,  on  M?\S(t).  Therefore,  the  aux¬ 
iliary  sources  (11)  may  only  differ  from  zero  in  the  transition  region  S(t)\S£(t). 
Note  also  that  the  overall  smoothness  of  the  original  solution  E(x,t),  H(x,t), 
as  well  as  that  of  the  multiplier  p(x,  t ),  guarantee  that  jm ,  j,  pm,  and  p  will  be 
smooth  compactly  supported  functions.  In  Figure  1,  we  schematically  depict 
the  geometry  of  the  region  on  which  the  auxiliary  sources  are  defined. 


Fig.  1.  Schematic  geometry  of  the  auxiliary  sources  region. 


Similarly  to  the  actual  physical  currents  and  charges,  the  electric  auxiliary 
sources  j(x,  t )  and  p{x ,  t)  of  (11)  do  satisfy  the  continuity  equation  (6),  as  can 
be  verified  by  direct  substitution.  The  other  pair  of  auxiliary  sources  jm  and 
pm  of  (11)  can  formally  be  interpreted  as  magnetic  currents  and  charges.  They 
are,  however,  nonexistent  in  nature  and  as  such  shall  be  considered  artifacts 
with  no  actual  physical  meaning.  They  are  only  introduced  for  intermediate 
derivation  purposes,  and  will  never  be  a  part  of  any  final  result.  The  magnetic 
sources  jm(x,t )  and  pm(x,t )  also  satisfy  the  continuity  equation  (6). 

The  auxiliary  problem  can  now  be  formulated:  To  integrate  the  Maxwell  equa¬ 
tions  driven  by  the  source  terms  (11)  on  the  entire  space  subject  to  the  homo¬ 
geneous  initial  conditions;  all  equations  are  inhomogeneous  rather  than  only 
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the  second  pair,  as  in  system  (4).  The  auxiliary  problem  is  linear  everywhere 
on  M3.  Under  the  assumption  of  sufficient  regularity,  its  solution  is  unique. 
The  regularity  is  guaranteed  by  the  requirement  that  the  transition  from  the 
exterior  solution  on  R3\S(t )  to  zero  on  Se(t)  be  smooth,  see  (9)  and  (10),  so 
that  the  auxiliary  sources  (11)  are  smooth  as  well.  Uniqueness  for  the  auxil¬ 
iary  problem  implies  that  its  solution  will  coincide  on  the  entire  space  with  the 
modified  fields  E{x,t)  and  H(x,t)  of  (10).  This,  in  particular,  means  that  on 
M?\S(t)  the  auxiliary  solution  will  coincide  with  the  original  exterior  solution 
E  (x,t),  H(x,t).  Note  that  we  do  not  need  to  know  the  latter  on  M.3\S(t)  in 
order  to  obtain  the  source  terms  jm(x,t),  j(x,t),  pm(x,t),  and  p(x,t)  of  (11) 
that  drive  the  auxiliary  problem.  We  only  need  to  be  sure  that  it  satisfies  the 
homogeneous  Maxwell  equations  (8)  and  that  transition  (10)  is  smooth. 

Altogether,  we  have  effectively  split  the  original  problem  into  two:  The  lin¬ 
ear  auxiliary  problem  that  needs  to  be  solved  on  the  entire  space,  and  the 
interior  problem  on  S(t)  that  will  be  integrated  with  the  external  boundary 
data  provided  by  the  solution  of  the  auxiliary  problem.  In  their  own  turn,  the 
sources  that  drive  the  auxiliary  problem  depend  on  the  interior  solution.  The 
algorithm  for  solving  the  auxiliary  problem  will  rely  on  the  presence  of  lacunae 
in  its  solutions,  see  Sections  2  and  3.  This  algorithm  will  allow  one  to  keep 
the  size  of  the  auxiliary  domain  finite,  computational  expenses  per  unit  time 
interval  fixed,  and  accuracy  non- deteriorating  over  long  time  intervals.  The 
lacunae-based  algorithm  is  described  in  Section  6.  Its  implementation  will,  in 
fact,  require  introducing  certain  modifications  to  the  formulation  of  the  aux¬ 
iliary  problem.  These  rather  minor  modifications  will  affect  the  construction 
of  the  auxiliary  sources  so  that  to  guarantee  existence  of  the  lacunae. 


5  The  Transverse-Magnetic  Subsystem 


From  here  on,  we  employ  cylindrical  coordinates  (r,  6 ,  z)  in  space  and  assume 
axial  symmetry  so  that  all  the  unknown  quantities  will  depend  only  on  r,  z, 
and  the  time  t,  and  J^(  •  )  =  0.  This  facilitates  the  split  of  system  (4)  into 
two  independent  subsystems  such  that  each  one  governs  only  three  out  of  the 
total  of  six  held  components.  The  subsystem  that  connects  Eg,  Hr,  and  Hz  is 
referred  to  as  transverse-magnetic  (TM);  it  contains  three  unsteady  equations 

1  dEe  _  fdHr  _  d Hz\ 
c  dt  \  dz  dr  ) 

1  dHr  _  dEg 
c  dt  dz 
1  dHz  |  ld(rE0) 
c  dt  r  dr 


47T  . 

— Je 
c 


0 

0 


(13a) 
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supplemented  by  the  steady-state  equation 


1  d(rHr)  |  OIL 
r  dr  dz 


(13b) 


The  latter  obviously  corresponds  to  V  •  H  —  0.  Note  that  the  second  steady- 
state  equation  of  (4),  V  •  E  —  An p,  does  not  appear  in  the  TM  system  because 
=  0.  It  is  rather  a  part  of  the  other  subsystem  known  as  transverse- 
electric  (TE);  the  latter  will  not  be  considered  in  the  current  paper. 


It  is  also  known  that  in  electrodynamics  the  six  field  components  are  not  com¬ 
pletely  independent,  they  are  determined  by  the  vector  and  scalar  potentials  A 
and  (j),  respectively.  The  potentials,  in  their  own  turn,  are  not  defined  uniquely. 
A  particular  gauge  [23]  that  we  are  choosing  here  to  specify  the  potentials  is 
H  =  V  x  A,  E  =  —  0  =  0,  which  translates  into 

dAg  Id  ( rAe )  1  dAe 

’  Hz  ~ - b - ’  ~ - AT 

dz  r  dr  c  dt 

In  other  words,  everything  is,  in  fact,  controlled  by  one  quantity  Ag,  which  is 
the  angular  component  of  the  vector  potential. 


It  is  instrumental  to  compare  equations  (13a),  (13b),  and  (14)  with  the  cylin- 
drically  symmetric  linearized  Euler  equations  that  appear  in  acoustics.  If  we 
introduce  the  new  variables:  rAg  :=  <p,  —rHr  :=  w,  rHz  :=  u,  rEg  :=  - p,  then 
system  (13a)  becomes 


1  dp 
c 2  dt 


1  d[ru) 
r  dr 

du 

dt 

dw 

dt 


2 

— u 
r 


dw 

dz 

dp 

dr 

dp 

dz 


An 

- rje 

c 


0 

0 


and  equation  (13b)  transforms  into 


du  dw 
dz  dr 


(15a) 


(15b) 


Equality  (15b)  can  be  interpreted  as  the  condition  of  zero  vorticity,  assuming 
that  u  and  w  are  the  radial  and  axial  components,  respectively,  of  the  velocity 
vector.  In  that  case,  (15b)  obviously  implies  conservativeness  of  the  velocity 
field,  which  is  sufficient  for  existence  of  the  lacunae  in  the  solutions  of  the 
acoustics  system,  see  [3].  The  gauge  (14)  is  then  equivalent  to 


(16) 
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where  p  should  be  thought  of  as  the  acoustic  pressure,  and  ip  —  the  velocity 
potential.  Equalities  (15b)  and  (16)  do  look,  in  fact,  exactly  the  same  as  the 
corresponding  relations  in  acoustics.  Let  us  emphasize,  however,  that  in  spite 
of  the  substantial  similarities,  system  (15a)  is  not  quite  identical  to  the  actual 
acoustics  system.  The  difference  resides  in  the  boxed  term  in  the  first  equation 
of  (15a),  which  is  not  present  in  the  genuine  acoustics.  In  fact,  the  two  systems 
-  two-dimensional  acoustics  and  TM  Maxwell  —  look  completely  identical 
only  in  the  Cartesian  coordinates.  Still,  the  similarities  that  we  have  identified 
for  the  cylindrical  coordinates  will  help  us  build  the  numerical  algorithm  for 
the  Maxwell  equations  that  will  share  certain  components  with  the  previously 
built  procedure  for  the  acoustics  system,  see  [3]. 

Next,  following  the  discussion  of  Section  4,  we  assume  that  the  homogeneous 
version  of  equations  (13a),  (13b)  [cf.  formulae  (8)]  holds  outside  the  domain 
S(t).  Clearly,  in  the  cylindrically  symmetric  case  S(t)  can  be  defined  as  a 
two-dimensional  shape  of  variables  (r,z),  while  the  actual  three-dimensional 
domain  will  be  given  by  the  corresponding  body  of  revolution.  The  motion  of 
the  domain  must  obviously  be  aligned  with  the  z-axis  [cf.  formulae  (7)]: 

w0  =  w0 (t),  z0  =  Zo(t)  =  f  w0{r)dT,  x0  =  y0  =  0  (17) 

Jo 

Inside  the  bounded  region  S(t),  some  possibly  complex  processes  may  be  going 
on  that  manifest  themselves  by  the  radiation  of  electromagnetic  waves  in  the 
far  held.  We  assume  the  overall  unique  solvability,  introduce  the  multiplier 
//  =  n(r,z,t)  as  in  (9),  and  modify  the  solution  on  M3  according  to  (10): 

Eg(r,  z,  t )  i — >  n(r ,  z,  t)E0(r,  z,  t )  =  Eg(r,  z,  t ) 

Hr(r,  z,  t )  i — >  p(r ,  z,  t)Hr(r ,  z,  t)  =  Hr{r ,  z,  t)  (18) 

Hz(r,  z,  t )  i — >  z,  t)Hz(r,  z,  t )  =  Hz(r,  z,  t ) 


Then,  applying  the  differential  operators  of  (13a),  (13b)  to  the  modified  held 
components  of  (18)  we  obtain  the  auxiliary  RHSs  [cf.  formula  (11)]: 


47 r- 
Jmr 
C 

47 r- 
7  m  2 


471 -pm 


1  dEg  _  fdHr  _  OIL  \ 
c  dt  \  dz  dr  J 

1  d  Hr  _  dEg 
c  dt  dz 
1  dHz  1  d  ( rEg ) 
c  dt  r  dr 
1  d  ( rHr )  dHz 
r  dr  dz 


(19) 


that  may  differ  from  zero  only  on  S(t)\Se(t).  The  non-physical  magnetic  cur- 
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(20) 


rents  and  charges  always  satisfy  the  continuity  equation  [cf.  formula  (6)]: 

dpm  1  d(rj  m  r  l  +  dl  m  z  _  q 
dt  r  dr  dz 

obtained  by  straightforward  differentiation  of  (19).  We  note  that  in  the  TM 
context  there  is  no  continuity  equation  for  the  electric  currents  and  charges. 
The  reason  is  that  the  only  type  of  electric  sources  that  appear  in  the  TM 
mode  is  the  current  component  jg,  for  which  we  obviously  have  ^  =  0.  The 
electric  continuity  equation  is,  in  fact,  a  part  of  the  TE  subsystem. 


There  is  also  an  alternative  procedure  for  obtaining  the  auxiliary  source  terms, 
besides  the  one  given  by  formulae  (18),  (19).  Instead  of  applying  the  multi¬ 
plier  p(r,  z,t, )  directly  to  all  the  field  components,  as  done  in  (18),  let  us  first 
reconstruct  Aq  in  the  transition  region  S(t)\Se(t)  using  the  first  two  equa¬ 
tions  of  (14).  Then,  the  multiplier  is  applied  to  Eg  and  Ae,  and  the  modified 
magnetic  hied  is  obtained  by  differentiating  the  modified  potential: 


Ee(r,  z,  t )  i — >  p(r,  z,  t)Eg(r,  z,  t )  =  Eg(r,  z,  t ) 

tt  dAe  _ld(rA0) 

Hr  ,  Hz  on  S(tj\Se(tj 

dz  r  dr 

Ag(r ,  2:,  t)  i — *  n(r ,  z,  t)Ae{r ,  z,  t )  =  Ag(r,  z,  t ) 

dAg 


(21) 


Hr(r ,  z,  t) 


dz 


=  Hr  (r ,  z,  t) 


1  d  ( rAe )  . 

Hz(r,  z,  t )  i — > - \ =  Ez[r,  z,  t ) 


r  dr 


Replacing  the  straightforward  modification  procedure  (18)  by  the  new  scheme 
(21)  that  involves  reconstruction  of  Ag  induces  important  simplifications  in 
the  definition  of  the  auxiliary  RHSs  (19).  Namely,  by  adopting  (21)  we  will 

guarantee  that  pm  =  0,  i.e.,  +  =  0.  Besides,  the  continuity  equation 

(20)  will  degenerate  and  reduce  to  1 =  0. 


Returning  to  the  issue  of  comparison  between  electromagnetics  and  acoustics, 
we  note  that  in  three-dimensional  acoustics  there  is  no  such  a  thing  as  the 
continuity  equation  for  the  source  terms.  In  two  dimensions,  however,  and 
in  particular,  in  the  (r,  z)  coordinates,  an  equation  similar  to  (20)  can  be 
written  that  would  connect  the  vorticity,  see  formula  (15b),  with  the  curl  of 
the  forcing  that  drives  the  momentum  equation.  In  case  the  acoustic  velocity 
hied  is  conservative,  this  equation  degenerates  similarly  to  the  aforementioned 
degeneration  of  the  continuity  equation  for  the  magnetic  sources. 


We  should  also  mention  that  the  foregoing  construction  of  the  auxiliary  sources 
is  by  no  means  the  only  possible  one.  In  fact,  all  we  need  is  a  smooth  exten¬ 
sion  of  the  exterior  solution  inwards  that  would  transition  to  zero  at  a  given 
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“depth”  e,  and  would  also  maintain  certain  properties  described  in  the  follow¬ 
ing  Section  6.  There  may  be  different  ways  of  obtaining  this  extension,  not 
necessarily  based  on  applying  the  multiplier  //  to  the  interior  solution  on  the 
transition  region. 


6  Lacunae-Based  Integration 


As  formulated  in  Section  4,  the  auxiliary  problem  consists  of  solving  the 
Maxwell  equations  driven  by  the  source  terms  (19)  that  are  concentrated  in¬ 
side  S(t),  or  more  precisely,  on  S(t)\Se(t).  In  other  words,  we  need  to  integrate 
the  system  of  equations: 


Id  Ee 
c  dt 


dHr 

dHz\ 

47T  - 

dz 

dr  ) 

- Je 

c 

ldHr 

dEe  _ 

47T  - 

c  dt 

dz 

—  Jm  r 

C 

b  1  d(rEe)  _ 

47T  - 
Jra  z 

c  dt  r 
1  d  (■ rH , 
r  dr 


+ 


dz 


=  47T  pn 


(22) 


for  t  >  0  subject  to  the  homogeneous  initial  conditions.  The  source  terms  in 
(22)  may  be  defined  by  either  (18),  (19)  or  (21),  (19).  The  solution  will  need 
to  be  found  on  a  domain  slightly  larger  than  S(t)  for  all  t  >  0.  Given  <5  >  0, 
we  define  this  new  domain  as  Ss(t)  =  {(cc,  t)\  dist  (x,  S(t))  <  5};  obviously, 
S£(t)  C  S(t)  C  Sg(t),  see  Figure  1.  We  emphasize  that  the  sources  (19)  are  not 
independent.  They  can  only  be  obtained  in  the  context  of  problem  decomposi¬ 
tion,  as  described  in  Section  4.  In  the  current  section  however,  we  will  describe 
the  solution  methodology  for  the  auxiliary  problem  as  if  those  sources  have 
simply  been  given.  In  the  following  Section  7,  this  solution  methodology  will 
be  applied  in  the  decomposition  framework  for  actually  setting  the  ABCs. 


The  definition  of  the  source  terms  (18),  (19)  or  (21),  (19),  along  with  the 
unique  solvability,  imply  that  the  solution  of  system  (22)  will  coincide  with 
the  modified  fields  (18)  or  (21),  respectively.  This,  in  particular,  means  that 
outside  S(t)  it  will  coincide  with  the  non-modified  fields,  as  the  multiplier  /i 
is  identically  equal  to  one  there,  see  (9). 


Let  us  now  introduce  a  partition  of  unity  on  the  semi-infinite  interval  t  >  0: 

OO 

Vf  >  0  :  ^2  0(f  —  aTi)  —  1  (23) 

i= o 
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where  T  >  0  and  \  <  a  <  1  are  two  parameters,  and  0  =  ©(£)  is  a  smooth, 
even,  “hat” -type  function  with  supp0(t)  C  [— f,  f]: 


m 


|i. 

1  -0(<rT- t), 

©H), 


0  <  t  <  (<r-  |)T 
(a  -  |)T<  t  <  f 
t  <  0 


(24) 


Then,  we  can  partition  the  sources  of  (22)  and  obtain  for  i  =  0, 1,  2, . 


j0l)(r,  z ,  t)  =  je(r ,  2,  t)0(t  -  oTi), 
j$.(r >  A  *)  =  Jmr(r,  2,  t)0(t  -  o-Ti), 
Jm2(r’  A  *)  =  im2(r,  2,  t)0(t  -  oTi), 
P^(r,2,i)  =  pm(r,z,t)Q(t  -  oTi), 


supp jg\r,z,t)  C  g,: 
suppjW(r,2,i)  C  Q, 
suppj^  (r,2,t)  c  q?- 
supppW(r,2,i)  c  g, 


(25) 


where  according  to  (24)  each  domain  Qi  is  compact  in  both  space  and  time: 

Qi  =  {(r,z,t)\(r,z)  e  S(t),  (ai  -  ^)T  <t<(ai  +  ^)t}  (26) 


Next,  we  would  like  to  consider  a  collection  of  individual  subsystems  driven 
by  the  partitioned  RHSs  of  (25)  [cf.  formulae  (22)]: 


IdEf  _  f  8H -W  _  0if«\ 

c  dt  \  dz  dr  ) 

1  dH®  _  dEf 
c  dt  dz 

1  dHf>  1  d  ( rE 
c  <9t  r  dr 
1  d  ( rH  W)  3if(0 
r  dr  dz 


4:7Tf) 


W 

m 


(27a) 


For  each  i  =  0,1,2,...,  the  corresponding  system  (27a)  needs  to  be  solved  for 
t  >  (ai  —  |)T,  subject  to  the  homogeneous  initial  conditions. 


We  need  to  emphasize  though  that  the  split  (25),  (27a)  may  not  always  be 
legitimate.  Indeed,  as  can  be  generally  shown,  see  [23],  and  as  has  been  indi¬ 
cated  when  obtaining  the  source  terms  (19),  every  RHS  to  a  Maxwell  system 
must  always  satisfy  the  continuity  equation  of  type  (20).  It  is  obvious  that 
partitioning  (25)  may  disrupt  that,  because  multiplication  by  a  function  of 
time  Q(t  —  aTi)  will  affect  the  time  derivative  of  ph).  A  natural  remedy  would 
be  to  use  the  definition  of  the  auxiliary  sources  by  means  of  formulae  (21),  (19) 
that  involve  the  reconstruction  of  Ag  on  S(t)\Se(t),  rather  than  the  general 
definition  (18),  (19).  As  has  been  mentioned,  by  using  (21),  (19)  we  guarantee 
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that  pm  =  0  and  also  independently  V  •  jm  =  0,  so  that  the  continuity  equa¬ 
tion  (20)  is  satisfied  identically  both  before  the  split  (25)  and  after  the  split, 
i.e.,  individually  for  each  i  —  0, 1,  2, ... .  As  such,  we  will  always  be  assuming 
hereafter  that  formulae  (21),  (19)  have  been  used  for  defining  the  auxiliary 
sources,  so  that  system  (27a)  reduces  to 


1  d  Ef 

c  dt 


1 

c 


'dH® 

dH®\ 

dz 

dr  J 

T3e 

1  dH^ 

dEf 

4:77 '■■(i) 

c  dt 

dz 

: - 

c  Jmr 

'(*)  i  d 

(r£f )  _ 

47 

7 - 1 - 

t  r 

dr 

: - ?W 

c  Jmz 

1  d  ( rH b 


dr 


dHf 

dz 


=  0 


(27b) 


a 


where  V  •  =  -  - 


(rM 

\  1  J  mr 


+ 


dtil  - 


=  0. 


dr  ,  Qz  —  The  RHSs  of  each  system  (27b)  are 
compactly  supported  in  space-time  on  the  domain  Qt,  see  (26),  and  each  of 
those  systems  is  now  uniquely  solvable.  As  (23)  is  a  partition  of  unity,  we  have: 


OO 

jo(r ,  z,  t)  =  ^2j0%\r,z,t),  jmr(r ,  z,  t)  =  ^jW(r,z,t), 

i= 0  i=0 


oo 

Jm  Z(r,z,t)  =  ^j£(r,^,t) 
i= 0 


and  because  of  the  linear  superposition,  the  solution  Eg(r,z,t),  Hr(r,z,t), 
Hz(r,z,t)  of  system  (22)  with  the  source  terms  defined  by  (21),  (19)  and 
subject  to  the  homogeneous  initial  conditions  at  t  —  0,  can  be  expanded  in 
terms  of  the  solutions  to  systems  (27b): 


Eg(r ,  z,  t)  =  ^2  Eg\r,  z,  t),  Hr(r,z,t )  = 


i= 0 


i= 0 


H2(r ,  z,  t)  =  Y^HP(r,z,t) 


(28) 


i=0 


The  series  (28)  are  formally  infinite.  However,  for  any  t  >  0  and  (r,z)  G  S(t) 
each  will,  in  fact,  contain  only  a  finite  fixed  number  of  nonzero  terms.  Indeed, 
due  to  the  causality  for  a  given  t  >  0  and  all  (cri  —  |)T  >  t,  i.e.,  all  i  > 
(|;  +  |)/cr,  we  will  have  Egl\r ,  z,  t )  =  0,  H^(r,  z,  t )  =  0,  and  H^(r,  z,  t)  —  0  on 
the  entire  space  M3.  Moreover,  since  the  RHSs  of  system  (27b)  are  compactly 
supported  in  space  and  time  on  the  domain  Qt  of  (26),  solutions  of  this  system 
are  going  to  have  lacunae,  as  shown  in  Section  3.  Let  us  now  denote  the 
moment  of  inception  of  source  ff  i  as  t^  =  ( ai  —  \)T  and  the  moment  of  its 

cessation  as  t\l>  =  (cri  +  ^)T.  Recall  also  that  diarn  S(t)  =  d,  which  means  that 
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cliam  S$(t)  =  d  +  25,  and  that  the  maximum  speed  of  motion  of  the  domain 
S(t)  is  k  <  c.  Then,  formula  (3)  implies  that  no  later  than 


_4.(i)  i  ,  d  +  25  +  (t[  —  tf)(c  +  k) 
“  0  +  +  c-k 

=  (ai  -  hr  +  d  +  2,5 ± _r(e  +  k)  = 

2  c  k 


+  Tint 


(29) 


all  of  the  domain  Ss(t)  will  fall  into  the  lacuna  of  system  (27b)  (see  [1,2]  for 
a  more  detailed  argument),  and  will  remain  inside  the  lacuna  continuously 
thereafter,  i.e.,  for  all  t  >  tf.  In  other  words,  for  the  solution  of  system 
(27b),  we  will  have  Ef(r,  z,  t)  =  0,  Hf(r,  z,  t )  =  0,  and  Hf(r,  z,  t)  —  0  when 
(r,z)  E  Ss(t)  and  t  >  t%\  where  is  defined  by  formula  (29).  Alternatively, 
this  means  that  for  any  t  >  0  and  all  i  <  (*~Tnt+|)/cr,  where  Tint  =  d+2S+T{<-+k) 

as  in  formula  (29),  the  terms  Eg\r,  z,t),  H^\r,z,t)  and  H^(r,z,t)  in  the 
series  (28)  will  be  equal  to  zero  for  (r,z)  E  Sg(t).  Consequently,  for  all  t  >  0 
and  (r,  z)  E  S$(t)  we  can  replace  expansions  (28)  with 


P  2 


Eo(r,  z,t)=  Y,  Ee\r,  z ,  t)%  Hr(r,  z,t)  =  Y  Hv\ri  zi  *), 


P  2 


l=P  1 


l=P  1 


P  2 


(30) 


l=P  1 


where  p1  =  +  \)/cr 


P  2  = 


(T  +  l)/cr  +1,  and  [  •  ]  denotes  the  integer 


part.  This  implies  that  for  any  t  >  0  and  (r,z)  E  S$(t )  the  number  of  terms 
p  =  P2  —  Pi  +  1  in  the  sum  (30),  and  as  such,  the  number  of  non-zero  terms  in 
the  expansion  (28),  may  not  exceed  [^W|  +  2.  Most  important,  this  number  p 


does  not  increase  as  the  time  t  elapses,  because  the  interval  Tint  introduced  in 
(29)  depends  only  on  the  partition  size  T  for  the  sources,  the  diameter  of  the 
domain,  the  propagation  speed,  and  the  maximum  speed  of  motion. 


Next,  we  realize  that  for  any  given  i  =  0,1,2,...  no  wave  can  travel  in  space 
further  away  than  the  distance  cTint  from  the  boundary  of  the  domain  S(t q')) 
during  the  time  interval  Tmt.  This  means  that  we  will  also  have  Eg  \r,z,t )  = 
0,  H^(r,z,t)  =  0,  and  H^\r,z,t)  =  0  for  dist  [(r,  z),  Sit^)]  >  cTint  and 


t. 


W 


o  <  t  <  4^-  As  such,  instead  of  the  free  unobstructed  space  outside  S(t) 


we  may  consider  outer  boundaries  with  arbitrary  (reflecting)  properties  for 
solving  each  of  the  problems  (27b).  As  long  as  none  of  these  boundaries  is 
located  closer  than  cTint  to  S(t^),  the  solution  of  (27b)  inside  Ss(t)  is  not 
going  to  feel  their  presence  for  4°  <  t  <  tf. 


In  fact,  instead  of  requiring  that  no  wave  may  reach  an  outer  boundary  before 
t  =  tf  we  can  introduce  a  weaker  requirement  that  no  reflected  wave  may 
come  back  and  reach  Ss(t)  before  t  =  tf.  The  latter  consideration  easily 
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translates  into  the  following  estimates  for  the  minimal  distances  between  the 
domain  S(t^)  and  the  allowed  location  of  any  reflecting  boundary,  [1,2]: 


7 


c  +  k 


T, 


int  5 


JD  _  _rpm 

2  int 


(31) 


Note  that  the  minimal  distances  in  the  directions  z  and  r  are  not  the  same,  see 
(31),  because  the  motion  of  the  domain  with  the  maximum  speed  k  is  allowed 
in  the  direction  z  only,  see  (17). 


Altogether,  the  solution  of  each  system  (27b)  driven  by  the  sources  compactly 
supported  on  the  domain  Qt  of  (26)  and  subject  to  the  homogeneous  initial 
conditions  at  t  —  t$  \  can  be  obtained  on  Sg(t)  for  all  t  >  as  follows.  First, 
this  solution  can  be  formally  considered  equal  to  zero  for  0  <  t  <  t$.  Next, 
on  the  time  interval  <  t  <  t%\  see  formula  (29),  system  (27b)  can  be 
integrated  on  the  rectangular  auxiliary  domain  with  the  sides  Z  and  R: 

Z  —  d  +  28  +  2Zmin  =  d  +  25  +  (c  +  k)Rnt 

R  =  d  +  8  +  R.min  =  d  +  5  +  -Tint 

centered  around  S(t^).  This  is  going  to  yield  the  correct  solution  inside  Sg(t). 
Finally,  for  all  t  >  the  solution  on  Sg(t)  will  be  equal  to  zero  again  because 
all  the  waves  will  have  left  the  domain  by  t  —  (lacuna). 


Let  us  now  return  to  the  representation  of  the  solution  of  system  (22)  with 
the  RHSs  given  by  (21),  (19)  in  the  form  of  finite  sums  (30).  Assume  that  the 
i-th  component  of  the  solution  is  obtained  by  integrating  the  corresponding 
system  (27b)  on  the  appropriate  domain  of  size  (32)  by  means  of  a  consistent 
and  stable  finite-difference  scheme.  System  (27b)  only  needs  to  be  integrated 
from  =  ( ai  —  |)T  till  =  t ®  +  Tint  because  for  all  subsequent  moments 
of  time  its  solution  on  Sg(t)  will  be  equal  to  zero.  Consequently,  the  following 
convergence  estimates  will  hold  for  to  *  <  t  <  1 2  *  and  (r,  z)  £  Ss(t): 

||st)(r,2,t)  -S<'‘-i)(r,2,«)||  <K,h‘ 

)  -  H^(r,z,t)\\  <  K,h‘  (33) 

where  a  is  the  order  of  convergence,  h  denotes  the  generic  grid  size,  and 
the  functions  E(,h’l\r,z,t ),  z,t),  and  H^(r,z,t)  denote  the  discrete 

solution  of  system  (27b)  for  a  given  i.  The  stability  constant  Kt  on  the  right- 
hand  side  of  each  inequality  (33)  does  not  depend  on  h,  but  may  depend  on 
3o\  j£.,  and  jg,  as  well  as  on  Tint. 

We  emphasize  that  the  quantity  Tint  does  not  depend  on  i.  Moreover,  it  is 
natural  to  assume  that  the  derivatives  of  the  functions  jg\r,  z,  t ),  j^r(jr:  z,  t), 
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and  j${r,z,t)  are  uniformly  bounded  with  respect  to  i.  In  this  case,  there 
will  be  a  universal,  i.e. ,  i-independent,  constant  K  =  K(jg,  jmr,  jmr,Tmt )  such 
that  Vi  =  0, 1,  2, . . .  :  K,  <  K.  Then,  using  representations  (30)  one  can 
easily  transform  the  individual  convergence  estimates  (33)  into  the  overall 
temporally  uniform  grid  convergence  estimate  for  Eg,  Hr,  and  Hz  that  would 
hold  for  t  >  0  and  (r,z)  G  Ss(t): 

II Eg(r,z,t)  -  E[eh\r,z,t)\\  <  pKha 

|| Hr{r,z,t)  -  Hlh\r,z,t)\\  <  pKha  (34) 

|| Hz(r,  z,  t )  -  H^h)(r,  z,  t)||  <  pKha 

In  formulae  (34),  E{eh)  =  £  ifW  =  £  H^\  and  ifW  =  £  H^\ 

i=pi  i=pi  i=pi 

cf.  formulae  (30).  A  detailed  proof  of  the  temporally  uniform  grid  convergence 
for  the  wave  equation  can  be  found  in  our  previous  work  [1], 

In  practical  terms,  the  temporally  uniform  grid  convergence  (34)  means  that 
the  accuracy  of  the  numerical  solution  of  system  (22)  with  the  RHSs  (21),  (19), 
will  not  deteriorate  even  on  arbitrarily  long  time  intervals,  if  the  integration 
is  performed  using  lacunae,  i.e.,  by  solving  the  set  of  systems  (27b)  and  then 
employing  representation  (30).  In  other  words,  one  should  expect  that  there  will 
be  no  long-time  error  buildup.  This  is,  in  fact,  a  key  distinction  between  the 
foregoing  lacunae-based  algorithm  and  traditional  time-marching  techniques 
that  may  be  applied  to  computing  the  unsteady  wave  fields. 

Indeed,  the  phenomenon  of  error  accumulation  during  long  runs  is  well-known 
in  the  context  of  computational  methods  for  time-dependent  problems.  At  the 
analysis  stage,  it  manifests  itself  by  the  increase  of  the  stability  constants  with 
time,  which  is,  in  fact,  equivalent  to  non-uniformity  of  the  grid  convergence. 
All  conventional  discrete  approximations  that  can  be  and  are  used  in  modern 
numerical  methods  are  known  to  suffer  from  this  deficiency. 

As  we  have  seen,  the  lacunae-based  algorithm  allows  us  to  circumvent  this 
difficulty,  because  even  though  for  any  individual  system  (27b)  the  stability 
constant  K,  will  still  depend  on  the  integration  interval,  the  duration  of  this 
interval  Tint  is  fixed  and  non-increasing  for  all  i  =  0, 1,2, ... .  Moreover,  the 
number  of  components,  see  (30),  needed  to  obtain  the  solution  on  Ss(t)  at  any 
given  moment  of  time  t  >  0  is  also  fixed  and  non-increasing,  which  altogether 
translates  into  the  temporally  uniform  convergence  (34). 

We  also  note  that  if  system  (22)  with  the  RHSs  given  by  (21),  (19)  were 
to  be  integrated  on  some  large  interval  [0,  7finai]  using  a  straightforward  time- 
marching  algorithm,  it  would  have  also  required  a  large  domain  in  space,  of  the 
size  roughly  2cTflna].  This  is  typically  not  feasible.  On  the  other  hand,  imple¬ 
mentation  of  the  lacunae-based  algorithm  allows  us  to  perform  the  integration 
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on  the  auxiliary  domain  of  a  fixed  and  non-increasing  size  determined  by  for¬ 
mulae  (32).  This  is  a  key  consideration  that  allows  us  to  use  the  lacunae-based 
algorithm  for  setting  the  ABCs,  see  Section  7,  i.e.,  for  obtaining  computational 
closures  when  the  original  unbounded  domain  of  the  problem  is  truncated. 

It  is  important  to  mention  that  smoothness  plays  a  key  role  in  the  design 
of  the  lacunae-based  algorithm.  In  particular,  the  function  0(f)  of  (24)  that 
helps  us  build  the  partition  of  unity  (23)  has  to  be  chosen  sufficiently  smooth 
so  that  the  dependence  of  the  stability  constants  K,  on  the  properties  of 
individual  RHSs  jjf\  ,  and  ,  see  (33),  be  not  “worse”  than  that  in  the 
original  scheme  with  non-partitioned  source  terms.  In  this  paper,  we  leave  out 
the  detailed  analysis  that  involves  the  quantitative  smoothness  characteristics, 
and  instead  refer  the  reader  to  our  previous  work  [1], 

Let  us  now  address  some  implementation  issues  for  the  lacunae-based  algo¬ 
rithm.  In  theory,  system  (27b)  for  every  given  i  =  0, 1,  2, . . .  should  be  inte¬ 
grated  on  its  own  auxiliary  region  of  the  size  given  by  (32),  centered  around 
the  domain  S(t^),  where  t^  =  (ai  —  \)T  and  the  location  of  S(t^)  is,  in 
turn,  determined  by  the  reference  point  r  =  0,  z  =  z0(t o'*),  see  formula  (17). 
It  is,  however,  more  convenient  to  consider  periodic  boundary  conditions  with 
the  period  Z  in  the  coordinate  direction  z.  This  essentially  means  that  the 
motion  (17)  should  be  interpreted  as  motion  on  the  circle.  In  so  doing,  all  sys¬ 
tems  (27b)  can  basically  be  solved  on  one  and  the  same  domain  of  variables 
(r,z):  [0,  A]  x  [—Z/2,Z/2],  Indeed,  it  does  not  matter  where  on  the  period 
the  “initial”  domain  S(t^)  is  located  for  every  i  —  0, 1,  2, ... .  The  boundary 
conditions  in  the  r  direction  are  discussed  in  Section  8. 

We  now  recast  the  discrete  version  of  formulae  (30)  in  the  form  of  a  difference: 

V  2  P2  pi-1  _ 

Eeh\r,z,t )  =  £  Egh,l\r,  z,  t)  =  £^M(r,z,f)  -  ^  E^l\r,z,t) 

i=pi  i= 0  i=0 

P2  P2  Pl-1 

Hih\r,z,t)  =  £  =  £ffM(r„M)  -  £  H<rh-'\r,z,t)  (35) 

i=pi  i= 0  i=0 

P2  P2  Pl-1 

-£«?■»(«•,*,»)  -  £  H^(r,z,t) 

i=pi  i= 0  i=0 

Besides  their  equivalence  to  (30),  formulae  (35)  can  be  given  the  following 

pi-i 

useful  interpretation.  The  subtracted  quantities  57  •  •  •  hr  (35)  are  assumed 

i=0 

zero  inside  the  lacuna.  More  precisely,  those  are  small  quantities  that  con¬ 
verge  to  zero  when  the  grid  is  refined.  Therefore,  their  subtraction  does  not 
change  the  solution  on  Ss(t),  it  rather  helps  halt  the  error  buildup  that  would 
continue  otherwise,  i.e.,  if  the  computation  of  these  terms  were  to  go  further 
on.  However,  outside  S$(t)  the  subtracted  terms  correspond  to  the  waves  that 
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have  left  the  domain  of  interest.  In  the  presence  of  reflecting  outer  boundaries, 

Pi-i 

the  subtraction  of  X]  ■  •  ■  in  (35)  helps  prevent  those  waves  from  coming  back 
i= o 

and  re-entering  the  domain  of  interest  Sg(t)  and  as  such,  “contaminating”  the 
solution  there. 

Existence  of  the  upper  limit  i  =  p-2  in  the  summation  (35)  is  due  to  the 
causality  which  is  always  a  factor,  and  has  nothing  to  do  with  the  lacunae. 

P2 

Therefore,  each  minuend  X  •  •  •  in  formulae  (35)  could  have  simply  been  ob- 

i= o 

tained  by  a  straightforward  time-marching  of  system  (22)  on  the  interval  [0,  t), 
with  absolutely  no  regard  to  either  the  partition  (25)  or  split  systems  (27b). 
Of  course,  neither  of  the  full  quantities  Egh\r,  z,  t),  (r,  z,  t ),  or  Hih)(r,z,t ) 

can  be  obtained  by  only  marching.  To  properly  account  for  the  presence  of 
pi— i 

the  subtrahends  X  '  ' '  in  formulae  (35),  let  us  first  symbolically  write  down 
i= o 

the  time- marching  scheme  that  would  apply  to  (22)  and  (27b): 

EP(r,z,  t  +  At)  =  £(£?'(-,  t),  H<h>(;  t),  fff  (-,*), M;  *)) 
H?\r,z,t  +  At)  =  Hr{EP(-tt),Hlh\;t);Jmr(;t))  (36) 

We  assume  that  the  boundary  conditions  are  built  into  the  operators  S,  TLr 
and  Ttz.  Note,  scheme  (36)  is  chosen  two-level  explicit  for  simplicity  only,  this 
is  by  no  means  a  limitation,  and  the  analysis  for  multi-level  schemes  can  be 
found  in  [1],  Consider  now  a  particular  moment  of  time  t  that  corresponds 
to  the  change  in  the  lower  summation  limit  in  formulae  (35)  from  i  =  p±  to 
i  —  pi  +  1,  i.e.,  such  t  that 


Pi  T  1  Pi 


Combining  formulae  (35)  and  (36),  we  will  then  have  for  this  particular  t: 

EP(;t  +  At)  =  £(£#■>(.,«),  t ),  H±%,t) ,M;  t)) 

-  E,pn>(;,t  +  At) 

Hih\;t  +  At)  =  nr(E(eh\;t),  HP);  t),jmr( ;  t))  -  t  +  At)  W 

HP(;t  +  At)=nz{EP(;t),HP(;t),jmz(;t))-Hil‘’-'\;t  +  At) 

In  other  words,  when  the  current  moment  of  time  t  satisfies  the  “switch”  condi¬ 
tion  (37),  the  terms  Egh’Pl\r,  z,t  +  At),  E[(h’Pl\r,  z,t  +  At),  and  i/|h’Pl)(r,  z,  t+ 
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At)  need  to  be  explicitly  subtracted  from  the  respective  overall  expressions, 
see  (38),  on  top  of  the  standard  time-marching  step  as  per  (36).  This  basi¬ 
cally  amounts  to  the  required  change  of  the  upper  summation  limit  in  both 
subtrahends  of  formulae  (35)  from  p\  —  1  to  p\.  We  will  also  assume  here¬ 
after  that  similarly  to  the  original  differential  equations  (22),  the  scheme  (36) 
satisfies  the  linear  superposition  principle.  Then,  the  next  time  step  after  the 
one  defined  by  formulae  (38),  i.e.,  the  step  t  +  At  \ — >  t  +  2 At,  shall  only  be 
done  by  marching  (36).  Indeed,  the  subtracted  quantities  Egh,pi\r,  z,t  +  At), 
h/Wpi)(r,  z,t  +  At),  and  H^h’Pl\r,  z,t  +  At)  will  carry  over  to  all  the  steps  that 
follow  (38)  due  to  the  linearity.  The  genuine  “unperturbed”  marching  can  thus 
continue  till  the  next  switching  moment  t,  i.e.,  till  condition  (37)  is  satisfied 
by  pi  +  2  and  p\  +  1.  At  this  moment,  the  quantities  Egl'Pl+1\r,  z,t  +  At), 
H(h’pi+1\r,z,t  +  At),  and  H^h’pi+1\r,  z,t  +  At)  will  need  to  be  subtracted, 
and  then  the  procedure  will  cyclically  repeat  itself. 

We  therefore  conclude  that  the  lacunae-based  algorithm  can  be  implemented 
as  a  conventional  time-marching  procedure  supplemented  by  repeated  sub¬ 
traction  of  the  retarded  terms,  see  (38).  The  subtraction  moments  can  be 
determined  ahead  of  time  and  are  separated  from  one  another  by  equal  time 
increments.  The  subtracted  terms  are  legitimately  called  retarded  because 
for  a  given  moment  of  time  t  that  satisfies  (37),  they  are  generated  by  the 
RHSs  that  are  active  in  the  past,  on  the  time  interval  [7qPi\  +  T]  = 
[{(Jpi  —  2 ) T,  (crpi  +  |)T]  =  [t— Tint,  t  —  Tint+T].  Of  course,  the  actual  subtracted 
quantities  Egl,pi+1\r,  z,  t+At),  H^h,Pl+1\r,  z,  t+At),  and  H^h’Pl+1\r,  z,  t+At) 
shall  be  re-computed  for  every  p\  independently  of  the  primary  time-marching 
procedure.  This  is  done  by  means  of  the  same  scheme  (36)  applied  to  the  cor¬ 
responding  system  (27b)  on  the  interval  [t  —  Tint,t]- 


7  Setting  the  ABCs 


Assume  that  there  is  a  space-time  grid  NxT,  on  which  a  discrete  approxima¬ 
tion  to  the  auxiliary  problem  is  built.  The  spatial  grid  N  consists  of  the  nodes 
n ,  whereas  the  temporal  grid  T  is  composed  of  the  time  levels  /  =  0, 1,  2, ... . 
Note  that  in  the  case  of  a  staggered  scheme  (see  the  example  in  Section  8) 
each  “formal”  node  (n,  l)  may,  in  fact,  correspond  to  a  collection  of  several 
space-time  locations,  each  associated  with  a  different  variable,  that  are  shifted 
with  respect  to  one  another  in  a  particular  way.  The  grid  N  is  actually  intro¬ 
duced  on  the  auxiliary  domain  [0,  A]  x  \—Z/ 2,  Z/ 2]  of  size  (32).  We  emphasize 
that  the  grid  N  does  not  have  to  offer  any  special  geometric  features  to  ac¬ 
commodate  the  shape  of  the  domain  S(t).  It  is  most  convenient  to  simply  use 
a  uniform  Cartesian  grid.  We  also  note  that  the  original  problem  solved  inside 
S(t)  does  not  have  to  be  approximated  on  the  same  grid.  In  the  most  gen- 
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eral  situation,  we  will  have  different  grids  for  the  interior  problem  and  for  the 
exterior /auxiliary  problem.  Then,  in  the  transition  region  S(t)\S£(t )  we  may 
need  to  employ  a  chimera-type  grid  strategy,  i.e. ,  interpolate  in-between  the 
overlapping  grids.  For  the  analysis  in  the  current  paper,  however,  we  will  sim¬ 
ply  assume  that  the  quantities  from  the  interior  solution  are  already  available 
on  the  grid  sub-domain  {Ff  x  T}  n  {S(t)\S£(t)}. 

Let  us  denote  by  Ffz,  l  =  0,1,2,...,  the  corresponding  time  levels  of  the 

grid  N  x  T,  and  by  §/  the  stencil  of  the  scheme  associated  with  the  node 

(n,l)  G  N  x  T.  For  simplicity,  we  will  assume  that  the  auxiliary  problem  (22), 

(21),  (19)  is  time-marched  by  an  explicit  scheme,  and  that  the  node  (n,  /)  G  E>ln 

is  the  actual  upper-level  node  on  the  m  +  1-level  stencil,  i.e.,  E>ln  fl  {Ff  x  T}  C 

(h/  U  F/_1  U  . . .  U  F/_m}  and  S^ON1  =  (n,l).  In  Section  6,  we  have  assumed 

m  =  1,  see  formula  (36).  Denote  by  F/+  the  sub-grid  of  F/  that  belongs  to  the 

interior  domain,  i.e.,  Ff^  —  NlD  S(tl),  where  tl  is  the  Z-th  time  level  in  actual 

units  (“seconds”);  in  the  simplest  case  tl  =  l  At.  Then,  introduce  the  sum  of 

the  interior  sub-grids  for  all  time  levels:  Ff+  =  N+  U  N+  U  Ff+  U  . . .  C  N  x  T. 

Finally,  consider  a  somewhat  larger  sub-grid  of  Ff  x  T:  Ff+ =  U  S*  which 

(n,/)eN+ 

is  simply  a  composition  of  all  the  stencils  Sln  obtained  when  the  upper-level 
node  (n,l)  sweeps  the  grid  Ff+;  obviously,  Ff+  C  Ff+.  The  part  of  the  grid 
Ff+  that  does  not  belong  to  Ff+  is  called  the  grid  boundary  and  is  denoted 
7  =  Ff+\Ff+.  We  will  require  that  the  domain  Ss(tl)  be  chosen  so  that  on 
every  time  level  tl,  l  —  0, 1,  2, ... ,  all  of  the  grid  N+  belong  to  this  domain: 
F/+  C  Ss(t l);  equivalently,  we  may  require  that  C  Ss(tl).  We  note  that  the 
grid  boundary  7  is  a  narrow  fringe  of  grid  nodes  that  follows  the  geometry  of 
dS(t).  Therefore,  the  parameter  6  may  be  chosen  small,  on  the  order  of  a  few 
grid  sizes  depending  on  the  specific  structure  of  the  stencil  E>ln.  We  also  note 
that  our  construction  of  the  grid  boundary  7  is  a  simplified  interpretation  of 
the  rigorous  general  definition  that  is  used  for  obtaining  the  discrete  Calderon’s 
potentials  and  boundary  projection  operators  [24];  the  latter  are  used  in  [16] 
as  a  universal  apparatus  for  setting  the  ABCs. 

We  will  now  describe  one  time  step  of  the  combined  time-marching  algorithm 
that  involves  setting  the  lacunae-based  ABCs.  In  so  doing,  we  will  assume  that 
all  time  steps  are  identical  and  as  such,  will  provide  an  inductive  description 
of  the  algorithm.  In  other  words,  all  the  operations  performed  on  a  given  time 
step  will  be  assumed  to  have  been  performed  on  all  the  previous  steps  as  well. 

Suppose  that  we  have  obtained  the  solution  for  up  to  a  given  time  level  l. 
This  means  that  the  solution  is  known  not  only  on  the  interior  domain,  but 
also  on  the  grid  boundary  —  at  the  levels  J1,  ■  ■  ■  ,  7*_m  that  are  immediately 
needed  for  advancing  the  next  time  step,  as  well  as  at  all  the  preceding  levels. 
In  particular,  one  may  think  about  starting  the  computation  from  the  known 
(homogeneous)  initial  conditions.  First,  we  make  one  interior  time  step  and 
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obtain  the  discrete  solution  everywhere  inside  S(t)  including  the  transition 
region,  i.e.,  the  grid  area  Nl+1  D{S  (tl+1)\Se(tl+1)} .  Then,  we  perform  the  mod¬ 
ification  (21)  in  the  discrete  framework,  which  involves,  among  other  things, 
reconstruction  of  Aq  on  the  grid  in  the  transition  region  S(tl+1)\S£(tl+1).  The 
most  straightforward  approach  to  do  that  is  contour  integration  along  the  grid 
lines;  it  has  proven  quite  robust  in  our  simulations,  see  Section  8,  taking  in 
to  account,  of  course,  that  the  potential  Ae  only  needs  to  be  reconstructed  on 
a  narrow  near-boundary  strip.  Having  gotten  the  modified  quantities  (21)  on 
the  grid  sub-domain  Nl+1  fl  {S  (tl+1)\S£(tl+1)} ,  we  apply  the  discrete  version 
of  (19)  [introduced  in  Section  8]  and  obtain  one  more  time  level  of  the  discrete 
auxiliary  RHSs.  If  the  scheme  written  on  the  stencil  approximates  system 
(22)  with  the  design  accuracy  at  some  node  (n  —  tiq,1  —  l0)  G  §(,  (clearly, 
l0  <  m),  then  the  discrete  RHSs  should  obviously  be  referred  to  this  same 
node  (n  —  no,  l  —  l0 )  and  as  such,  advancing  the  interior  solution  till  (/  +  1) 
would  mean  building  the  auxiliary  RHSs  till  (/  +  1  —  /o).  Next,  we  make  one 
time  step  for  the  auxiliary  problem  with  these  newly  updated  RHSs,  and  ob¬ 
tain  its  solution  on  Ni+1.  Since  we  have  chosen  5  >  0  so  that  C  Sg(tl+1), 
we  determine  that  the  solution  to  the  auxiliary  problem  will,  in  particular,  be 
available  on  ^l+l .  This  concludes  one  full  time  step,  because  once  we  know  the 
solution  on  all  time  levels  up  to  (/  +  1)  everywhere  including  the  grid  boundary, 
we  can  advance  the  next  interior  time  step,  etc. 

We  should  emphasize  that  the  auxiliary  problem  is  solved  by  the  laeunae- 
based  methodology  of  Section  6.  The  latter  includes  cyclic  subtractions  (38) 
of  the  retarded  terms  defined  by  the  partition  (25)  on  top  of  the  straightfor¬ 
ward  time-marching.  It  is  important  to  realize  that  once  a  particular  compo¬ 
nent  has  been  subtracted,  there  will  never  be  a  need  to  analyze/incorporate 
it  again  in  the  course  of  computation.  In  other  words,  the  subtracted  terms 
are  completely  disregarded  from  the  moment  of  subtraction  further  on.  There¬ 
fore,  the  corresponding  partition  elements  (25)  of  the  auxiliary  RHSs  can  be 
disregarded  as  well.  Consequently,  even  when  integrating  over  arbitrarily  long 
time  intervals,  we  only  need  to  keep  a  finite  amount  of  the  past  information, 
namely,  the  auxiliary  RHSs  (19)  defined  on  the  interval  of  duration  Tint,  see 
(29),  that  immediately  precedes  the  current  moment  of  time.  This  makes  the 
extent  of  temporal  nonlocality  of  the  proposed  ABCs  fixed  and  limited. 

Summarizing  on  the  implementation  strategy,  we  time-march  the  combined 
problem  by  performing  the  alternating  interior /auxiliary  steps.  To  advance 
an  interior  step,  we  need  the  data  on  the  grid  boundary  7,  i.e.,  at  those  nodes 
of  the  scheme  stencil  that  do  not  belong  to  the  interior  domain.  Providing 
these  missing  data  basically  means  setting  the  discrete  ABCs.  This  is  done 
in  practice  by  computing  the  solution  to  the  auxiliary  problem  on  the  finite 
domain  [0,  R]  x  [-Z/ 2,  Z/ 2]  of  size  (32),  which,  in  particular,  means  obtaining 
it  at  the  grid  boundary  7  as  well.  The  latter  solution  is,  in  turn,  driven  by  the 
source  terms  that  are  generated  based  on  the  interior  solution. 
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In  the  end  of  the  section,  we  outline  the  important  design  features  of  the  new 
unsteady  ABCs.  The  boundary  conditions  are  built  directly  in  the  discrete 
framework  and  can  be  constructed  for  any  given  scheme.  The  computational 
expenses  per  unit  time  still  remain  fixed  and  non-growing.  The  outer  boundary 
becomes  completely  transparent  for  all  the  outgoing  waves,  and  the  accuracy 
of  the  boundary  treatment  is  fully  determined  by  that  of  the  scheme  em¬ 
ployed.  The  fixed  and  non-growing  expenses  per  time  step  are  accounted  for 
by  using  the  lacunae-based  integration.  Besides,  the  lacunae-based  approach 
guarantees  that  the  ABCs,  while  being  nonlocal,  still  have  only  limited  and 
non-growing  extent  of  temporal  nonlocality.  In  addition  to  that,  the  accuracy 
of  the  boundary  treatment  per  se  will  not  deteriorate  as  the  time  elapses. 
Finally,  the  ABCs  have  the  capability  of  handling  moving  domains,  includ¬ 
ing  those  engaged  in  an  accelerated  motion.  In  the  following  Section  8,  we 
corroborate  the  foregoing  design  properties  of  the  proposed  ABCs  by  actual 
numerical  examples. 


8  Numerical  Demonstrations 


For  the  purpose  of  constructing  a  discretization,  it  will  be  necessary  to  write 
the  TM  Maxwell  system  on  the  axis  of  symmetry  r  —  0.  First  of  all,  we  note 
that  all  the  unknown  quantities  must  be  continuous  and  bounded  everywhere. 
For  the  vector  components  Eg  and  Hr  this  consideration,  along  with  the  axial 
symmetry,  imply  that  both  Eg  and  Hr  must  be  equal  to  zero  at  r  =  0.  The 
same  will  obviously  be  true  for  Ag.  For  the  vector  component  Hz ,  which  is 
parallel  to  the  axis  of  symmetry,  the  same  arguments  lead  to  the  conclusion 
that  =  0.  Next,  using  Taylor’s  expansion  for  r  <C  1,  we  have:  Eg(r ,  •)  = 

ar  r*=0 


E'g( 0,  -)r  +  o(r)  and  consequently,  = 


that 


1  9(rEe) 
r  dr 


r=0 


O  dEp 

dr 


.  Therefore,  for  r 

r= 0 


l  d(r2  E'JO,-))  .  /-■  \  ,-i 

— -  +  o(l),  which  means 

=  0  system  (22)  transforms  into: 


Ee(0,z,t)  =0 
Hr( 0,  z,t)  =  0 


1  dHz 

c  dt 

2  dJk 

dr 


+  2 
■  + 


d  Eg 
dr 
dHz 
dz 


0 

0 


(39) 


which,  in  particular,  implies  that  the  RHSs  (21),  (19)  are  equal  to  zero  on  the 
axis.  The  gauge  definition  (14)  for  r  =  0  becomes: 


Hz  = 

dr 


(40) 
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The  domain  S(t)  is  chosen  as  a  ball  of  fixed  diameter  d  centered  on  the  z 
axis,  with  its  center  moving  according  to  (17).  Obviously,  as  the  motion  (17) 
is  aligned  with  the  direction  r  =  0,  it  does  not  break  the  axial  symmetry.  The 
value  for  the  diameter  that  we  take  is  d  —  1.8;  the  functions  wo(t),  and  zo(t) 
of  (17)  will  be  specified  later.  The  value  for  the  speed  of  light  is  chosen  c  =  1. 

The  auxiliary  domain  is  a  rectangle  [0,  R]  x  [— Z/2,  Z/2]  of  variables  (r,  z ),  with 
the  actual  sizes  R  =  t:  and  Z  =  2n.  We  make  sure  that  the  other  parameters 
a,  T,  and  5,  see  (29),  are  chosen  so  that  these  sizes  be  greater  or  equal  to  the 
corresponding  minimal  values  prescribed  by  (32).  The  boundary  conditions  in 
the  z  direction  are  periodic: 

Ee(r,z±  Z,t)  =E0(r,z,t ) 

Hr(r,  z  ±  Z,  t)  =  Hr(r,  z,  t )  (41a) 

Hz(r,z±  Z,t)  =Hz(r,z,t ) 


The  boundary  conditions  in  the  radial  direction  obviously  cannot  be  periodic 
because  of  the  geometry/symmetry  considerations.  At  r  =  0,  there  is  no  need 
for  special  boundary  treatment  at  all;  instead,  we  have  the  first  two  equations 

=  0.  At  r  =  R,  the  boundary  conditions  are  needed 

r=0 

for  the  homogeneous  counterpart  of  system  (22),  and  we  set: 


of  (39)  along  with  ^ 


d(rEe 


dr 


r=R 


=  0, 


d(rHr 


dr 


r=R 


=  0,  H~(R ,  z,t)  —  0 


(41b) 


We  note  that  the  particular  reflecting  properties  of  boundary  conditions  (41a) 
and  (41b)  are  basically  immaterial  from  the  standpoint  of  reconstructing  the 
infinite-domain  solution  on  S(t),  as  long  as  the  geometric  restrictions  discussed 
in  Section  6  are  honored.  We  only  need  to  make  sure  that  the  auxiliary  prob¬ 
lem  is  uniquely  solvable,  which  boundary  conditions  (41a),  (41b)  do  provide 
for.  We  leave  out  a  detailed  justification  of  the  latter  statement,  only  mention 
that  the  homogeneous  Robin  boundary  condition  (41b)  for  Eg  follows  from 
the  homogeneous  Dirichlet  boundary  condition  (41b)  for  Hz  combined  with 
the  homogeneous  counterpart  of  the  third  equation  of  (22);  and  the  homoge¬ 
neous  Robin  boundary  condition  (41b)  for  H r  follows  from  the  homogeneous 
Dirichlet  boundary  condition  (41b)  for  Hz  combined  with  the  homogeneous 
counterpart  of  the  fourth  (steady-state)  equation  of  (22). 


To  validate  the  proposed  numerical  method,  a  test  solution  is  needed.  Con¬ 
struction  of  the  test  solution  for  electromagnetics  will  be  somewhat  more  elab¬ 
orate  than  that  for  acoustics,  see  [3].  In  acoustics,  we  have  been  able  to  build 
the  test  solution  by  first  obtaining  the  velocity  potential,  which  is  a  genuine 
scalar  function.  In  the  context  of  electrodynamics,  the  quantity  Ag  that  de¬ 
termines  all  the  fields  according  to  (14)  is  rather  a  component  of  a  vector.  In 
the  TM  mode,  this  component  may  not  depend  on  the  polar  angle  9. 
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We  start  with  constructing  two  auxiliary  dipole  solutions  ipx  and  ifjy  of  the 
three-dimensional  scalar  wave  equation: 

l_d  _  A{V>X,  V’y}  =  x(t){Sx,Sy}(x,y,  z  -  z0(t))  (42) 

where  in  each  case  the  RHS  is  a  point  dipole  oriented  either  along  x  or  along  y, 
which  travels  according  to  the  law  (17)  and  has  a  time-dependent  amplitude 
x(t).  Solutions  of  equation  (42)  can  be  obtained  by  computing  convolutions 
with  the  fundamental  solution  of  the  wave  operator,  which  is  a  single  layer  of 
unit  magnitude  on  the  expanding  sphere  of  radius  ct,  see  [3]  and  [25,  Chapter 
7]  for  detail.  To  write  down  i[}x  and  in  a  convenient  form,  let  us  introduce  a 

new  function  f  =  r(t)  =f  \Jr2  +  (z  —  z0(t))2  and  denote  v  =  r(r)  +  ct,  where 
t  is  a  parameter.  Then, 


fp{x,y}(r,z,t)  = 


+ 


J- _ c{x,y}x'(T) 

4vr  [cr(r)  -  (z  -  20(t))w0(t)]2 

c{x,y}x(r) 

4:71 


lS=Ct 


_r(r)[cr(r)  -  (z  -  z0(t))w0(t )]2 

.o(j-^|)TO|T)  +  w2(t)  _(z_  Zo(t))<{t  - 


(43) 


cr  (t 


(z  -  ^0(r))w0(^)]3 


v=ct 


The  right-hand  side  of  equality  (43)  has  to  be  evaluated  for  v  =  ct,  i.e.,  for 

\jr2  +  (z  —  Zq{t))2  +  ct  =  ct  (44) 

Solution  r  of  the  algebraic  equation  (44)  gives  the  actual  retarded  moment  of 
time,  at  which  the  trajectory  of  the  traveling  point  source  intersects  the  lower 
portion  of  the  characteristic  cone  with  the  vertex  (r,z,t).  For  the  case  of  a 
uniform  motion,  equation  (44)  is  quadratic  with  respect  to  r  and  can  be  solved 
in  the  closed  form.  Its  root  r  <  t,  once  substituted  in  (43),  yields  solution  of 
the  wave  equation  that  can  also  be  obtained  via  the  Lorentz  transform.  This 
approach  was  adopted  in  [1,  2]  for  the  case  of  monopole  sources.  In  general 
however,  one  should  not  expect  to  be  able  to  solve  the  nonlinear  equation  (44) 
explicitly.  For  the  simulations  in  the  current  paper,  equation  (44)  is  solved 
numerically  using  the  Newton  method. 

Next,  we  know  that  the  vector  potential  A  satisfies  the  (vector)  wave  equation, 
see  [23].  Therefore,  its  Cartesian  components  should  individually  satisfy  the 
standard  scalar  wave  equations.  Let  us  now  define  Ax  =  ^  and  Ay  =  ^x, 
see  (42).  Then,  expressing  the  angular  component  of  the  potential  as  Aq  = 
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Ax  sin  9  +  Ay  cos  9  we  obtain: 


A0(r,z,t)  = 


cr 


4ir  [cr (r)  -  (z  -  2o(t))wo(t)]2  X  ^  + 

gfggg  +  »g(r)  -  B(2  -  2.(t)K(t) 

cr(r)  -  (2  -  Zo(t))wo(t) 


(45) 


v=ct 


Expression  (45)  obviously  gives  an  appropriate,  i.e.,  axially  symmetric,  A$.  To 
actually  evaluate  it  for  a  given  ( r,z,t ),  the  corresponding  retarded  moment  r 
needs  to  be  determined.  Assuming  that  k\  <  w0(t )  <  k2,  where  k\  and  k2  are 
known,  we  solve  equation  (44)  by  Newton’s  method  with  the  initial  guess  To 
given  by  the  solution  of  the  quadratic  equation:  4  r2  +  (z  —  l(ki  +  /c2)ro)2  = 
c(t-T0)2  that  satisfies  To  <  t.  For  the  numerical  demonstrations  in  this  paper, 
we  specify  the  same  law  of  motion  (17)  as  we  used  in  acoustics,  see  [3]: 


Zo(t) 


3 

+  T-S 


where  [  •  ]  denotes  the  integer  part,  as  before,  and  {  •  }  denotes  the  fractional 
part.  The  motion  (46)  consists  of  repeated  acceleration/deceleration  cycles 
of  duration  10,  so  that  during  each  cycle  the  source  travels  a  total  distance 
of  1  along  the  z  axis.  Both  the  velocity  w0 (t)  =  z'0(t )  and  the  acceleration 
w'0  (t)  =  Zo(t )  determined  by  (46)  are  continuous  functions  of  time,  and  0  = 
k\  <  wait)  <  k2  =  0.1875,  which  means  that  the  condition  max*  |w0(t)|  <  c  is 
met  because  c  —  1. 


Solution  (45)  for  the  potential  Ae  is  obviously  singular,  and  cannot  be  used 
directly  to  derive  the  quantities  needed  for  testing  the  numerical  procedure. 
To  remove  the  singularity,  we  introduce  a  smooth  function  G  =  G(r),  r  =  f(t ), 
such  that  G(r)  =  1  for  r  >  'y,  where  n  <  1,  and  also  such  that  d  =  0 
for  k  —  0, . . .  ,  4.  Then,  the  new  function  Ag(r,z,t )  •  G(r)  is  continuous  and 
bounded  everywhere.  By  differentiating  it,  we  define  the  components  of  the 
reference  solution  to  the  cylindrical  TM  system  [cf.  formulae  (14)]: 

Eo(r,z,t )  =  -  ^^A0(r,z,t)  ■  G(r)) 

Hr(r,z,t )  =  -  ^(A0(r,z,t)  ■  G(f ))  (47) 

Hz(r ,  z,  t)  =-^-(r-Ae(r,z,t)  ■  G(f)) 

On  the  axis  r  =  0,  we  use  the  first  two  relations  of  (39)  and  formula  (40).  Note 
that  Aq  of  (45)  is  a  retarded  potential,  and  therefore  differentiating  it  as  in  (47) 

4  Obtained  from  (44)  for  the  case  of  a  uniform  motion  with  the  speed  (Aq  +  k2)/ 2. 
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will  involve  implicit  differentiation  of  r  via  equation  (44).  In  contradistinction 
to  that,  the  definition  of  G(r)  does  not  involve  any  retardation. 

The  variables  Eg,  Hr,  and  Hz  of  (47)  satisfy  the  homogeneous  Maxwell  sys¬ 
tem  for  r  >  4r,  i.e.,  outside  a  smaller  ball  concentric  with  S(t).  The  value 
of  k  that  we  took  for  our  simulations  was  k  =  0.8.  Inside  this  smaller  ball, 
i.e.,  for  r  <  Eg,  Hr ,  and  Hz  of  (47)  satisfy  equations  (13a),  (13b)  with 
the  current  jg  obtained  by  directly  substituting  the  held  components  (47)  into 
the  first  equation  of  the  system.  Note,  when  removing  the  singularity  of  the 
solution  we  have  required  that  sufficiently  many  derivatives  of  G  be  equal  to 
zero  at  r  =  0;  consequently,  the  resulting  jg  will  have  no  singularities  either. 
Altogether,  we  have  therefore  obtained  a  reference  electromagnetic  solution, 
which  is  regular  everywhere,  and  which  can  be  said  to  be  generated  by  the  elec¬ 
tric  currents  concentrated  inside  the  moving  domain  S(t),  more  precisely,  on  a 
smaller  concentric  ball  of  diameter  nd.  Outside  of  this  smaller  ball,  our  refer¬ 
ence  solution  is  basically  a  system  of  unsteady  electromagnetic  waves  radiated 
by  and  propagating  away  from  a  moving  source.  We  are  going  to  reconstruct 
this  solution  numerically  on  the  domain  S(t),  and  set  the  discrete  laeunae- 
based  ABCs  on  its  outer  boundary  dS(t )  according  to  the  methodology  of 
Sections  7  and  6.  The  latter  will  include  formulating  the  interconnected  inte¬ 
rior  and  auxiliary  problems,  building  the  near-boundary  auxiliary  RHSs,  and 
integrating  the  auxiliary  problem  using  the  lacunae-based  algorithm. 

The  Maxwell  system  (13a)  is  approximated  numerically  on  the  Cartesian  grid 
of  variables  (r,  z): 


Vi  =  iAr,  A r  =  R/Nr ,  i  —  0, 1, . . .  ,  Nr 


Zj  =  jAz,  Az  =  Z/2Nz,  j  =  0,  ±1, 


±N, 


(48) 


using  a  second-order  staggered  finite-difference  scheme: 
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(49a) 


Note  that  we  do  not  include  the  discretized  equation  (13b)  into  the  integration 
scheme  (49a).  From  the  structure  of  system  (13a)  it  is  easy  to  conclude  that  if 
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equation  (13b)  holds  for  t  —  0,  which  is  always  the  case  for  the  homogeneous 
initial  conditions,  it  will  also  hold  for  all  subsequent  t  >  0.  This  property  is,  in 
fact,  inherited  by  the  staggered  discretization  (49a),  as  can  be  verified  directly, 
assuming  that  the  discrete  derivatives  for  (13b)  are  done  in  a  natural  way 
consistent  with  (49a).  Note  also  that  the  auxiliary  system  (27b)  is  discretized 
using  the  same  scheme  (49a).  The  only  difference  is  that  the  second  and  third 
equations  may  be  inhomogeneous,  the  corresponding  discrete  RHSs  j$  and 
j$  are  referred  to  the  grid  locations  (■ i,j  +  |,  /)  and  (i  +  respectively. 

The  discrete  version  of  the  last  equation  of  (27b)  does  not  have  to  be  a  part 
of  the  integration  scheme  either.  Indeed,  as  indicated  in  Section  7  all  the 
auxiliary  RHSs  are  obtained  directly  on  the  discrete  level,  i.e.,  by  applying  the 
corresponding  difference  operators  to  the  modified  fields  on  the  grid,  similarly 
to  (19).  In  so  doing,  if  the  modification  (21)  is  used,  then  the  appropriate 
discrete  version  of  the  degenerate  continuity  equation  1  d('rJd™r -  +  (^r  =  0  will 
guarantee  that  once  the  discretized  last  equation  of  (27b)  holds  for  the  initial 
data,  it  will  hold  for  all  subsequent  moments  of  time  as  well. 


Scheme  (49a)  is,  in  fact,  a  version  of  the  well-known  Yee  scheme  that  was 
originally  introduced  for  solving  the  Cartesian  Maxwell  equations,  see  [26]. 
Discretization  (49a)  is  valid  for  i  >  0,  i.e.,  away  from  the  axis  of  symmetry. 
On  the  axis,  we  discretize  equations  (39),  which  yields: 
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Periodic  boundary  conditions  (41a)  are  recast  on  the  grid  as  follows: 


(49b) 
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(50a) 


whereas  boundary  conditions  (41b)  are  discretized  as: 


rN-E«Nr,i  -  r^E'e  i 
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(50b) 


A  staggered  second-order  scheme  similar  to  (49a),  (49b)  has  been  used  for 
conducting  the  numerical  experiments  with  lacunae-based  ABCs  in  acoustics, 
see  [3].  There  are,  however,  important  differences,  not  only  in  the  structure  of 
the  system,  see  (15a),  but  also  in  how  the  variables  are  assigned  to  the  integer 
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and  semi-integer  grid  nodes,  and  how  the  continuous  and  discrete  boundary 
conditions  in  the  radial  direction  are  set. 


Four  our  simulations  of  the  Maxwell  equations,  we  have  used  a  sequence 
of  three  successively  more  fine  square  cell  grids,  A r  =  A z,  of  dimensions 
Nr  x  2 Nz  =  64  x  128,  128  x  256,  and  256  x  512.  The  Courant  stability  con¬ 
straint  was  applied  when  selecting  the  time  step  for  the  explicit  scheme  (49a), 
(49b).  The  grid  boundary  7  was  built  in  accordance  with  the  definition  given  in 
Section  7,  taking  into  account  that  scheme  (49a),  (49b)  is  staggered.  The  lat¬ 
ter  implies  that  we  effectively  have  three  independent  stencils  that  are  shifted 
with  respect  to  one  another  for  updating  one  electric  field  component  and  two 
magnetic  field  components.  Of  course,  the  grid  boundary  is  constructed  con¬ 
currently  with  the  actual  time-marching.  The  parameter  5  that  is  needed  to 
accommodate  the  width  of  the  grid  boundary  7  was  taken  5  =  |  Ar.  The  func¬ 
tions  O (t),  n(x,t),  and  G(r)  introduced  in  order  to  guarantee  smoothness  at 
all  the  stages  of  the  derivation,  are  constructed  as  piece-wise  polynomials  with 
four  continuous  derivatives  everywhere.  In  so  doing,  the  multiplier  fi(x,t),  see 
(9),  is  also  built  as  only  a  function  of  r.  The  amplitude  y(f),  which  is  a  part 
of  the  definition  of  the  potential  Ag,  see  (45),  was  chosen  in  the  form  of  a 
harmonic  oscillation  with  the  frequency  three  times  that  of  the  motion  cycles 
(46).  The  width  £  of  the  transition  region  S(t)\S£(t),  see  Figure  1,  varied  to 
demonstrate  different  aspects  of  the  algorithm  performance,  see  below.  The 
parameter  cr,  which  is  a  part  of  the  definition  of  the  partition  of  unity  (23), 
(24)  was  chosen  a  —  |.  The  actual  temporal  thickness  T  of  each  partition  ele¬ 
ment  was  calculated  “backwards”  from  (29)  and  (32),  considering  that  the  do¬ 
main  size,  the  maximum  motion  speed,  see  (46),  and  the  period  Z  are  known. 
When  marching  the  auxiliary  problem,  retarded  components  of  the  solution 
are  subtracted  according  to  (38).  Each  subtracted  component  is  recomputed 
as  solution  to  the  corresponding  problem  (27b).  For  a  given  i,  this  problem  is 
inhomogeneous  on  the  interval  [4444  =  [(cri  —  \)T,  (ai  +  |)T],  and  then  it 

remains  homogeneous  on  the  interval  [4444  =  [(°^  +  \)T,  ( cri  —  |)T  +  T-mt\. 
Therefore,  we  first  explicitly  time-march  this  system  on  its  interval  of  inhomo¬ 
geneity.  Then,  we  do  the  FFT  of  the  solution  in  the  x  direction  and  expansion 
with  respect  to  the  corresponding  eigenfunctions  (evaluated  numerically)  in 
the  r  direction,  which  allows  us  to  advance  the  homogeneous  solution  further 
till  4^  by  simply  raising  the  resulting  amplifications  factors  to  the  correspond¬ 
ing  powers.  We  note  that  in  so  doing  scheme  (49a)  effectively  gets  split  into 
three  discrete  wave  equations  for  the  respective  unknown  quantities.  The  rea¬ 
son  is  that  because  of  the  cylindrical  geometry  different  transforms  along  r 
are  needed  for  different  variables.  This  is  basically  the  only  instance  at  which 
reduction  to  a  set  of  independent  wave  equations  is  employed;  and  it  is  only 
necessitated  by  a  particular  choice  of  the  coordinate  system,  which  allows  us 
to  take  advantage  of  the  three-dimensional  effects  in  a  two-dimensional  com¬ 
putational  setting. 
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Lacunae-based  ABCs  for  unsteady  Maxwell 


Grid  convergence  for  electric  field,  accelerated  source  motion 


Fig.  2.  Grid  convergence  study  with  lacunae-based  ABCs,  e  =  lOAr. 

In  Figure  2  we  present  the  computational  error  for  Eg  as  it  depends  on  time 
on  all  three  grids  that  we  have  employed.  The  total  integration  time  was  100^, 
i.e.,  one  hundred  times  the  interval  required  for  the  waves  to  cross  the  domain. 
The  error  was  evaluated  in  the  maximum  norm  on  the  domain  S(t).  We  see 
that  the  algorithm  provides  for  no  long-term  error  buildup,  and  also  that 
it  displays  the  design  second-order  grid  convergence.  Error  profiles  for  both 
magnetic  hied  components  Hr  and  Hz  look  practically  the  same  and  we  are 
not  showing  them  here. 

Similarly  to  our  analysis  in  acoustics  [3],  we  will  now  discuss  how  the  perfor¬ 
mance  of  the  lacunae-based  ABCs  is  affected  by  the  key  parameter  involved 
-  the  width  of  the  transition  region  S(t)\S£(t )  defined  Section  4.  This  width 
£  obviously  reflects  on  how  well  the  smooth  multiplier  function  fi(x,t)  is  re¬ 
solved  on  the  grid,  and  as  such,  how  smooth  the  auxiliary  RHSs  (19)  will 
effectively  be.  The  latter,  in  turn,  affect  the  quality  of  the  discrete  lacunae, 
i.e.,  how  sharp  the  aft  fronts  of  the  waves  really  are  in  the  discrete  framework. 
This  is  important  because  every  time  a  retarded  component  is  subtracted, 
see  (38),  we  assume  that  what  is  being  subtracted  inside  the  lacuna,  i.e.,  on 
the  domain  S$(t),  is  zero,  or  more  precisely,  a  small  quantity  that  converges 
to  zero  with  the  refinement  of  the  grid;  the  convergence,  however,  hinges  on 
the  smoothness  of  the  source  terms.  Besides  having  a  potential  effect  on  the 
error  behavior,  the  width  of  the  transition  region  also  determines  how  many 
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Lacunae-based  ABCs  for  unsteady  Maxwell 


Grid  convergence  for  electric  field,  accelerated  source  motion 


Fig.  3.  Grid  convergence  study  with  lacunae-based  ABCs,  e  =  8A r. 


grid  nodes  are  needed  to  support  the  auxiliary  RHSs.  We  remind  that  those 
RHSs  basically  control  the  extent  of  temporal  nonlocality  of  the  lacunae-based 
ABCs.  The  algorithm  requires  keeping  them  on  the  interval  of  length  Tint  that 
immediately  precedes  the  current  moment  of  time,  and  as  such,  the  more 
narrow  the  transition  region  is,  the  less  additional  storage  is  needed. 


The  error  curves  on  Figure  2  were  obtained  for  the  transition  region  S(t)\S£(t) 
being  relatively  wide,  on  the  order  of  ten  grid  cell  sizes.  In  so  doing,  the  actual 
geometric  width  of  S(t)\S£(t )  decreases  with  the  refinement  of  the  grid.  We 
are,  however,  interested  to  find  out  what  happens  when  the  number  of  cells  in 
the  transition  region  also  decreases.  In  Figures  3,  4,  5,  and  6  we  are  showing  the 
Eg  error  profiles  for  the  width  of  the  transition  region  being  e  =  8,  6,  4,  and  2 
grid  cell  sizes,  respectively. 

We  observe  that  with  the  decrease  of  £  the  error  deteriorates,  which  is  natural 
to  expect.  We  also  notice,  though,  that  the  deterioration  is  more  visible  on  the 
coarser  grids,  whereas  on  the  finest  grid  that  we  have  employed,  256  x  512,  it 
is  much  slower.  Qualitatively,  this  is  the  same  type  of  behavior  as  we  have  seen 
in  acoustics,  see  [3].  There  is  also  an  important  difference.  In  acoustics,  when 
£  decreases,  all  the  error  profiles  start  to  grow  more  or  less  monotonically 
with  time,  see  [3].  The  smaller  the  £  the  faster  this  growth  is,  although  on 
the  finest  grid  it  is  not  as  fast  as  on  the  coarser  grids.  In  contradistinction 
to  that,  in  the  current  electromagnetic  context  the  error  profile  on  the  finest 
grid  always  remains  flat  on  some  initial  time  interval,  after  which  it  starts 
to  deteriorate.  The  extent  of  this  initial  interval  decreases  with  the  decrease 
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Fig.  4.  Grid  convergence  study  with  lacunae-based  ABCs,  e  =  6Ar. 

Lacunae-based  ABCs  for  unsteady  Maxwell 


Grid  convergence  for  electric  field,  accelerated  source  motion 


Fig.  5.  Grid  convergence  study  with  lacunae-based  ABCs,  e  =  4Ar. 


of  e,  but  even  for  the  narrowest  transition  region  that  we  have  considered, 
£  =  2 A r,  it  is  still  quite  substantial,  about  30  units,  i.e.,  30  times  the  time 
needed  for  the  waves  to  cross  the  domain,  see  Figure  6.  The  presence  of  this 
initial  flat  portion  of  the  error  profile  is  beneficial  as  it  essentially  means  that 
the  lacunae-based  ABCs  can  be  used  for  a  certain  period  of  time  with  no 
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Fig.  6.  Grid  convergence  study  with  lacunae-based  ABCs,  e  =  2Ar. 

deterioration  of  performance  even  when  the  auxiliary  RHSs  involved  are  not 
really  smooth.  We  note  that  having  only  two  grid  sizes  across  S(t)\Se(t)  does 
imply  that  there  is  no  smoothing  in  the  transition  region  at  all.  We  are  rather 
having  a  sharp  truncation,  and  effectively  substituting  an  equivalent  of  the 
Heaviside  function  on  the  grid  for  the  multiplier  n(x,t). 

As  of  yet,  we  cannot  offer  a  rigorous  theoretical  explanation  of  why  the  algo¬ 
rithm  appears  more  sensitive  to  the  quality  of  the  discrete  lacunae  on  coarser 
grids  than  on  the  fine  grid.  To  some  extent  this  is  counterintuitive  because  a 
typical  instability  would  rather  manifest  itself  by  a  rapid  deterioration  of  the 
solution  when  the  grid  is  refined.  We  can  only  “speculate”  that  the  observed 
behavior  has  to  do  with  the  actual  magnitude  of  those  discrete  “tails”  behind 
the  aft  fronts  of  the  waves  that  are  due  to  the  “imperfections”  in  the  auxiliary 
sources,  and  that  apparently  are  still  smaller  on  fine  grids.  Altogether,  this 
phenomenon  is  certainly  advantageous  for  computations  with  lacunae-based 
ABCs,  because  fine  grids  are  needed  for  high  overall  accuracy  anyway,  and 
at  the  same  time  they  will  allow  to  maintain  high  accuracy  of  the  boundary 
treatment  for  longer  periods  of  time. 


9  Discussion 


We  have  constructed  and  tested  the  algorithm  for  setting  highly- accurate 
global  artificial  boundary  conditions  for  the  computation  of  time-dependent 
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electromagnetic  waves.  This  work  extends  our  previous  approach  that  applies 
to  the  scalar  wave  equation  [2],  and  is  also  similar  to  the  ABCs  that  we  have 
recently  developed  for  unsteady  acoustics  [3].  The  algorithm  is  based  on  the 
presence  of  lacunae  (aft  fronts  of  the  waves)  in  the  three-dimensional  wave- 
type  solutions.  The  ABCs  are  obtained  directly  for  the  discrete  formulation 
of  the  problem  and  can  complement  any  consistent  and  stable  finite-difference 
scheme.  In  so  doing,  neither  a  rational  approximation  of  non-reflecting  kernels, 
nor  discretization  of  the  continuous  boundary  conditions  is  required,  see  the 
review  [5] .  The  extent  of  temporal  nonlocality  of  the  new  ABCs  appears  fixed 
and  limited,  and  this  is  not  a  result  of  any  approximation  but  rather  a  di¬ 
rect  consequence  of  the  fundamental  properties  of  the  solution.  The  proposed 
ABCs  can  handle  artificial  boundaries  of  irregular  shape  on  regular  grids  with 
no  fitting/ adaptation  needed.  Besides,  they  possess  a  unique  capability  of  be¬ 
ing  able  to  handle  boundaries  of  moving  computational  domains,  including 
the  case  of  accelerated  motion. 

The  key  idea  of  using  lacunae  for  computations  is  straightforward:  One  should 
not  continue  the  computation  inside  the  lacuna  once  the  solution  has  become 
zero  there  due  to  “natural  causes.”  Otherwise  the  error  may  unwarrantably 
build  up  on  one  hand,  and  on  the  other  hand,  the  extent  of  the  required 
temporal  pre-history  of  the  solution  may  grow  un-justifiably  high. 

Of  course,  there  are  many  technical  issues  that  need  to  be  addressed  before 
this  idea  can  actually  be  implemented  in  practice.  Those  relate  primarily  to 
what  to  do  with  the  waves  outside,  rather  than  inside,  the  domain  of  interest. 
In  the  framework  of  the  previous  analysis,  we  simply  allow  them  to  propagate 
a  certain  distance  away  till  they  get  reflected,  and  then  set  up  the  auxiliary 
domain  so  that  the  reflected  waves  do  not  reach  the  domain  of  interest  before 
it  completely  falls  inside  the  lacuna.  This,  however,  is  by  no  means  the  only 
possible  option.  In  fact,  any  treatment  of  the  outgoing  waves  that  would  pre¬ 
vent  the  reflections  from  re-entering  the  domain  of  interest  before  it  falls  into 
the  lacuna  will  be  appropriate.  At  the  same  time,  introducing  an  alternative 
to  the  approach  described  above  may  be  beneficial  from  the  standpoint  of  the 
overall  computational  cost.  For  example,  a  treatment  of  the  sponge  layer  type 
that  slows  down  the  outgoing  waves,  see,  e.g.,  [27,28],  may  allow  to  reduce  the 
size  of  the  auxiliary  domain.  Alternatively,  the  lacunae-based  approach  can 
be  combined  with  a  PML-based  treatment  for  the  waves  outside  the  domain 
of  interest,  see,  e.g.,  the  survey  papers  [29,30]. 

The  scope  of  numerical  demonstrations  presented  in  the  current  paper  is  lim¬ 
ited  to  the  cylindrically  symmetric  case.  The  reason  for  that  is  two-fold.  On 
one  hand,  by  doing  so  we  can  take  advantage  of  the  fundamentally  important 
three-dimensional  effects  in  an  essentially  two-dimensional  computational  set¬ 
ting.  We  have,  in  fact,  followed  this  strategy  in  our  previous  work  on  acous¬ 
tics  [3]  and  on  the  scalar  wave  equation  [2],  On  the  other  hand,  cylindrically 
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symmetric  TM  mode  helps  us  guarantee  that  the  partitioned  auxiliary  sources 
will  satisfy  the  correct  continuity  equation,  see  the  beginning  of  Section  6. 
For  that  purpose  we  also  need  to  use  a  particular  method  of  obtaining  the 
auxiliary  sources  that  involves  reconstruction  of  the  potential  component  Aq. 
Reconstruction  of  Aq  in  the  context  of  electromagnetics  is  similar  to  recon¬ 
struction  of  the  velocity  potential  that  we  have  used  in  acoustics,  see  [3].  In 
electromagnetics,  however,  it  is  needed  to  ensure  the  solvability  the  auxiliary 
sub-problems  with  the  partitioned  RHSs  (27b),  whereas  in  acoustics  existence 
of  the  potential  implies  existence  of  the  lacunae  in  the  solutions  of  the  auxil¬ 
iary  problem.  We  should  note  that  in  the  full-fledged  electromagnetics  context 
that  involves  all  six  field  components  (as  opposed  to  the  TM  subsystem  that 
governs  only  three  components)  there  will  be  both  the  magnetic  continuity 
equation  and  the  electric  continuity  equation  for  the  RHSs  that  will  need  to 
be  enforced  once  the  sources  are  partitioned  in  time.  At  this  moment  of  time, 
the  question  of  how  it  may  be  achieved  for  the  electric  part  remains  open. 
To  some  extent,  this  presents  an  obstacle  for  constructing  the  lacunae-based 
ABCs  for  the  Maxwell  equations  in  the  most  general  case.  Of  course,  this  ob¬ 
stacle  can  be  easily  overcome  by  simply  reducing  the  system  to  a  collection 
of  six  individual  wave  equations  for  the  field  components,  and  treating  each 
equation  independently  from  the  others  according  to  the  methodology  of  [2], 
This  always  remains  a  “fail-safe”  option;  it  may,  however,  be  more  expensive 
than  obtaining  the  ABCs  in  the  framework  of  one  first-order  system.  As  such, 
implementation  for  the  full  three-dimensional  Maxwell  system  will  be  in  the 
focus  of  the  future  study. 

For  the  current  cylindrically  symmetric  case,  numerical  experiments  conducted 
in  the  paper  fully  corroborate  the  theoretical  design  properties  of  the  method. 
Besides,  when  a  key  parameter  of  the  algorithm,  the  width  of  the  transition 
region,  changes  so  that  to  reduce  the  overall  memory  requirements,  the  per¬ 
formance  of  the  ABCs  deteriorates  in  an  expected  way.  However,  on  fine  grids 
this  deterioration  is  considerably  slower  than  that  on  coarse  grids.  Therefore, 
for  those  practical  problems,  for  which  typical  integration  times  are  sufficiently 
long  but  not  excessively  long,  this  deterioration  my  go  almost  unnoticed. 
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